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SECTION  I 
INTRODUCTION 


The  signing  of  the  Limited  Nuclear  Test  Ban  Treaty  in  1963  made  necessary  the 
development  of  methods  for  simulating  nuclear  airblast  and  groundshock 
effects.  The  High  Explosive  Simulation  Technique  (HEST)  was  developed  during 
the  mid-1960s  to  simulate  airblast  pressure  and  airblast-induced  groundshock 
loadings  on  near-surface  buried  structures  (Refs.  1  and  2).  HEST,  as 
originally  devised,  did  not  simulate  the  dynamic  airblast  effects  (reflected 
pressures,  etc.)  which  would  be  experienced  by  above-ground  structural 
targets.  The  Dynamic  Airblast  Simulator  (DABS)  was  developed  during  the 
mid-1970s  to  provide  for  the  testing  of  above-ground  missile  shelter  concepts 
(Ref.  3). 


The  development  of  these  and  other  nuclear  airblast  and  groundshock  simula¬ 
tion  methods  over  the  past  15  years  was  accomplished  by  primarily  experimental 
methods.  Empirically  derived  expressions  based  on  previous  experimental  data 
were  used  for  the  design  and  sizing  of  future  tests  (Ref.  4). 


As  the  requirements  for  more  accurate  simulation  of  a  wider  range  of  nuclear 
blast  environments  continued  to  grow,  these  experimental  and  empirical  design 
methods  became  less  adequate.  Certain  problems  continued  to  plague  each  of 


1.  Auld,  H.  E.,  D'Arcy,  G.  P.,  and  Leigh,  G.  G.,  Simulation  of  Airblast- 
Induoed  Ground  Motions  (Phase  I),  AFWL-TR-65-1 1 ,  Air  Force  Weapons  Labora¬ 
tory,  Kirtland  Air  Force  Base,  New  Mexico,  April  1965. 

2.  Auld,  H.  E.,  D'Arcy,  G.  P.,  and  Leigh,  G.  G.,  Simulation  of  Airblast- 
Indueed  Ground  Motions  (Phase  II),  AFWL-TR-65-26,  Vol .  I,  Air  Force  Weapons 
Laboratory,  Kirtland  Air  Force  Base,  New  Mexico,  April  1965. 

3.  Martens,  Daniel  P.,  and  Bradshaw,  Joel  C.,  Dynamic  Airblast  Simulator 
Parametria  Test  Series ,  Events  I-A ,  I-B,  I-C}  I-D ,  and  I-E  Data  Report, 

AFWL-TR-76-018,  Air  Force  Weapons  Laboratory,  Kirtland  Air  Force  Base, 

New  Mexico,  November  1976. 

4.  Bratton,  J.  L.,  and  Pratt,  H.  R.,  Simulation  of  Airblast- Induced  Ground 

Motions  (Phzse  HA),  AFWL-TR-66-85,  Air  Force  Weapons  Laboratory,  Kirt¬ 
land  Air  Force  Base,  Mew  Mexico,  October  1967. 
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these  simulators.  Not  only  did  variation  of  the  experimentally  derived 
design  parameters  fail  to  provide  solutions,  but  certain  phenomena  observed 
in  the  data  defied  explanation  by  purely  empirical  means. 


Over  the  past  15  years  substantial  progress  has  been  made  in  the  area 
of  first-principle  hydrodynamic  calculations.  Several  large  hydrodynamic 
computer  codes  (hydrocodes)  have  evolved  with  which  the  simulation  methods 
can  be  examined  in  greater  detail.  A  program  of  hydrocode  calculations  of 
the  HEST  and  DABS  simulators  was  undertaken  to  help  understand  the  basic 
explosion  mechanics  of  these  simulators  and  to  provide  improved  capabilities 
for  design  and  prediction. 

Some  difficulties  have  been  encountered  with  both  the  HEST  and  the  DABS 
simulators  which  could  not  be  resolved  through  experimental  or  empirical 
techniques.  In  HEST,  high  pressure  spikes  and  severe  oscillations  were 
often  encountered  near  the  beginning  of  the  pressure  pulse  (Fig.  1). 

Whether  these  were  caused  by  positioning  individual  strands  of  detonating 
cord  near  the  pressure  gages  or  by  larger  shock  oscillations  within  the  HEST 
explosion  cavity  is  not  known.  Tests  conducted  in  which  all  explosives  in 
the  cavity  were  initiated  nearly  simultaneously  resulted  in  greater  peak 
pressures  for  a  given  charge  density*  than  tests  in  which  the  detonation 
was  caused  to  propagate  through  the  explosion  cavity  at  some  specified  rate 
(sweeping  wave).  Finally  the  relation  between  the  charge  density  in  the  cav¬ 
ity  and  the  resultant  peak  pressure  was  not  well-defined  for  a  wide  range  of 
pressures  and  for  different  explosive  materials. 

As  in  HEST  the  DABS  technique  must  provide  waveforms  with  the  desired  peak 
pressure  and  impulse.  In  addition  DABS  must  generate  blast  waves  with  the 
appropriate  dynamic  pressure  characteristics.  These  conditions  alone  place 
severe  restrictions  on  the  design  of  DABS.  Coupled  with  the  requirement  for 
large-scale  testing  capability,  these  conditions  make  DABS  a  unique  challenge 
for  predictive  numerical  simulations. 


*  Charge  density--the  amount  of  explosive  material  per  unit  volume  of 
explosion  cavity. 
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HEST  Airblast  Pressure  Trace 


Difficulties  encountered  with  DABS  have  included  uncertainties  about  explosive 
charge  density  in  the  driver  section  versus  peak  pressures  and  impulse  at 
locations  down  range,  shock  interaction  from  ribs  in  the  trench  sections, 
shock  reflections  and  other  interactions  involving  the  contact  surface,* 
and  perturbation  of  the  flow  field  by  the  test  structure.  These  and  other 
design  considerations  led  to  the  need  for  hydrocode  calculations  to  support 
the  further  development  of  these  simulation  methods.  The  extent  to  which  the 
existing  hydrocodes  could  model  all  aspects  of  these  simulation  experiments 
was  not  known  nor  were  the  costs  of  performing  such  calculations.  Answers 
to  both  of  these  questions  were  to  be  obtained  during  this  technical  effort. 


♦Contact  surface--the  material  interface  between  the  detonation  products 
and  the  shocked  air. 
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SECTION  II 
CALCULATIONAL  TOOLS 


The  calculational  tools  used  during  this  technical  effort  included  the  one¬ 
dimensional  (1-D)  similarity  solution,  the  APOD  1  1/2-D  hydrocode,  the  SAP 
1-D  hydrocode  (both  Eulerian  and  Lagrangian  versions),  and  the  HULL  two- 
dimensional  (2-D)  hydrocode.  Each  is  briefly  described  below. 

SIMILARITY  SOLUTION 

The  1-D  similarity  solution  followed  the  method  of  Sakurai  (Ref.  5)  which 
provides  similarity  solutions  for  blast  waves  in  planar,  cylindrical,  and 
spherical  geometries.  A  small  code  was  written  to  computerize  this  method 
to  ease  the  solution  of  a  variety  of  problems.  Idealized  blast  wave  problems 
using  a  pure  energy  source  in  an  ideal  gas  can  be  calculated  out  to  the 
point  of  encounter  with  a  reflecting  surface.  The  similarity  solution, 
however,  will  not  handle  shock  reflections  and  thus  can  only  be  used  to 
estimate  the  very  early  portion  of  most  explosion  problems.  It  also  cannot 
be  used  to  calculate  the  detonation  of  high  explosive  materials  in  a  real 
gas  medium.  A  listing  (in  BASIC)  of  this  computer  code  is  provided  in 
Appendix  A. 


APOD  HYDROCODE 

APOD  is  a  1  1/2-D  Lagrangian  hydrodynamics  computer  code  developed  by  Mr.  Bob 
Port  of  Research  and  Development  Associates  (RDA),  Marina  Del  Rey, 

California.  The  finite  differencing  scheme  in  APOD  follows  the  method  of 
S.  K.  Godunov  (Ref.  6).  The  burn  option  is  somewhat  oversimplified  but  it 
does  allow  both  1 eft-to-right  and  right-to-left  burn  of  the  explosive  driver 
material.  No  mixing  of  material  is  allowed  by  this  code. 


5.  Sakurai,  A.,  "On  the  Propagation  and  Structure  of  the  Blast  Wave,  1," 

Journal  of  the  Phuoiaal  Jo.-’ct;-  of  Japan,  Vol .  8,  No.  5,  September- 
October  1  953. 

6.  Richtmyer,  Robert  D.  and  Morton,  K.  W.,  Piffcreme  Methods  for  Initial- 
Value  Fi'ohlcna,  2nd  Edition,  Interscience  Publishers. 
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A  real  air  equation  of  state  (EOS)  is  used  in  the  reaction  portion  of  the 
mesh.  However,  the  EOS  for  the  driver  section  is  a  modified  gamma-law  scheme 
which  approximates  the  driver  EOS  with  a  variable  gamma  dependent  upon  the 
ratio  of  change  in  cell  thickness.  Driver  mass  density  and  energy  density  are 
obtained  from  data  card  input  on  each  run. 

Cartesian  coordinates  are  employed  exclusively  in  APOD.  Reflective  and  trans¬ 
missive  options  for  the  right  boundary  condition  are  written  into  the  code, 
but  they  still  require  internal  code  modification  in  order  to  be  switched  from 
one  to  the  other.  The  left  boundary  condition,  or  backwall,  was  originally 
coded  as  an  "energy  sink,"  i.e.,  energy  was  allowed  to  pass  into  a  mesh  behind 
the  driver  section.  Because  thickness  of  the  cell  increased  with  distance 
behind  the  backwall  boundary,  pressure  waves  entering  this  region  were  never 
allowed  to  flow  toward  the  backwall.  In  a  recent  modification  of  APOD* 
the  backwall  is  now  treated  as  a  high-density  soil  whose  behavior  is  governed 
by  a  soil  lockup  model.** 

One  of  the  primary  advantages  of  APOD  is  the  scheme  for  the  expansion  of  the 
radial  wall.  This  scheme  leads  to  APOD's  being  referred  to  a  1  1/2-D  code. 

In  this  scheme  the  main  hydrocode  is  allowed  to  compute  the  flow  parameters 
along  the  main  axis.  The  pressure  within  a  given  cell  is  then  used  to 
compute  a  change  in  volume  of  that  cell  as  the  cross-sectional  area  expands 
due  to  the  pressure.  The  code  then  calculates  the  corresponding  change  in 
mass  density  and  loss  of  energy  due  to  expansion  of  the  cell  walls.  The 
same  soil  lockup  model  used  in  the  backwall  section  is  also  the  basis  of 
the  wal 1 -expansion  method. 

APOD  does  not  possess  a  restart  capability.  Fortunately,  because  it  is  a 
fast-running  code,  it  can  handle  all  but  extremely  long  calculations.  The 
plotting  options  for  APOD  are  limited.  There  is  no  provision  for  any 


*Chown,  W.,  and  Harrison,  B.  D. ,  Hydrocode  for  Calculating  Interaction  in 
Typical  DABS  Configurations ,  Task  Report,  UNM/CERF  AST-8,  Eric  H.  Wang  Civil 
Engineering  Research  Facility,  University  of  New  Mexico,  Albuquerque,  New 
Mexico,  February  1978. 

**Seusy,  F. ,  Lock-up  Iwpulsc  Code  Description ,  DED-A,  Technical  Memorandum, 
Air  Force  Weapons  Laboratory,  Kirtland  Air  Force  Base,  New  Mexico,  5 
November  1  976. 
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data  snapshots  or  profile  plots.  Histories  of  pressure  and  pressure-impulse 
times  are  overlaid  on  the  same  graph  for  a  maximum  of  10  stations. 

Organization  of  the  code  is  fixed  with  practically  no  flexibility  to  meet 
varying  cal cul ational  demands.  If  the  geometry  of  the  problem  changes,  major 
code  revisions  are  required.  Another  of  the  major  drawbacks  of  the  code  is 
the  reverse  sense  of  the  major  axis.  Originally  designed  for  a  specific 
purpose,  the  axis  runs  negative  to  the  right,  thus  reversing  the  sense  of  most 
of  the  logic,  a  factor  which  leads  to  considerable  confusion  when  internal 
code  tracing  is  attempted. 

SAP  HYDROCODE 

The  SAP  is  a  1-D  hydrocode  developed  at  the  Air  Force  Weapons  Laboratory  (AFUL) 
in  1966  by  Whitaker  et  al .  (Ref.  7).  In  1975  SAP  was  integrated  into  the  HULL 
system  which  used  the  SAIL  updating  management  system  (Ref.  8).  The  SAIL 
system  provides  numerous  variations  in  coding  by  the  appropriate  selection  of 
several  options.  This  allows  a  high  degree  of  versatility  in  specifying 
any  given  calculation. 

Two  finite  differencing  options  are  provided,  and  the  code  is  adaptable  with 
a  minimum  of  effort  to  the  differencing  scheme  desired.  The  first  method  is 
a  Lagrangian  method  described  in  detail  in  AFWL-TR-66-141  (Ref.  7).  The 
second  method  is  an  adaptation  of  the  HULL  differencing  scheme  which  offers 
either  Eulerian  or  Lagrangian  methods  and  is  described  further  in  AFWL-TR- 
76-183  (Ref.  9). 

7.  Whitaker.  W.  A.,  et  al..  Theoretical  Calculations  of  the  Phcnomenoloju  of 
HE  Detonations,  AFWL-TR-66-141,  Vol .  1,  Air  Force  Weapons  Laboratory, 
Kirtland  Air  Force  Base,  New  Mexico,  November  1966. 

8.  Graham,  D.  C.,  Gaby,  L.  P.,  II,  and  Rhodes,  C.  E.,  Jr.,  SAIL,  An  automated 

Approach  to  Cofi-jare  Development  and  Management ,  AFWL-TR-78-80,  Air  Force 
Weapons  Laboratory,  Kirtland  Air  Force  Base,  New  Mexico,  October  1976. 

9.  Fry,  M.  A.,  et  al  . ,  The  HULL  Hpdrodunariics  Computer  Code ,  AFWL-TR-76-1 83 , 
Air  Force  Weapons  Laboratory,  Kirtland  Air  Force  Base,  New  Mexico,  Sep¬ 
tember  1976. 


Cartesian,  cylindrical,  or  spherical  coordinate  systems  may  be  elected. 

SAP  may  be  run  on  smaller  calculations  with  or  without  tapes.  Because  SAP 
offers  a  restart  capability,  a  longer  calculation  may  be  continued  should  the 
calculation  be  interrupted  before  completion.  The  code  will  compute  flows  at 
any  angle  to  the  horizontal.  The  mesh  may  be  rezoned  periodically  so  that  the 
shock  is  always  located  interior  to  the  mesh. 

SAP  uses  a  real  atmosphere  that  is  stable  under  an  R"2  gravity  field  and  an 
air  equation  of  state  that  is  an  empirical  fit  to  Hilsenrath's  data  (Refs.  10 
and  11).  Six  additional  formulations  provide  EOSs  for  TNT,  PBX,  9404, 
Pentolite,  methane,  and  ammonium  nitrate  and  fuel  oil  with  provisions  for 
including  others  if  desired. 

A  continuous-burn  routine  based  on  Chapman-Jouget  theory  is  present  in  the 
code.  Any  explosive  material  may  be  used  if  EOSs  for  both  the  burned  and 
unburned  material  is  known.  This  routine  provides  the  conditions  existing 
in  the  gaseous  explosive  products  at  the  instant  the  detonation  wave  reaches 
the  surface  of  the  exploding  charge.  The  routines  work  for  forward  burn 
only.  In  other  words  the  process  must  be  initiated  by  starting  with  at  least 
one  cell  of  burned  material  at  the  left-most  perimeter  of  the  section  of  the 
mesh  containing  the  explosive  material. 

Some  other  features  of  SAP  are  as  follows.  First,  boundary  conditions  may  be 
specified  by  card  input  as  transmissive,  reflective,  or  a  specified  function 
of  time.  Second,  one  of  the  tapes  used  with  SAP  is  a  station  tape  which 
records  the  time  histories  of  hydrodynamic  variables  at  particular  (fixed) 
locations.  Next,  plotting  options  for  SAP  allow  graphs  of  pressure,  density, 
and  velocity  versus  radius  to  be  generated  at  specified  times  during  the 
calculation.  After  the  calculation  is  completed  a  station  plotter  provides 


10.  Hilsenrath,  J.,  Green,  M.  S.,  and  Beckett,  C.  W.,  ThewoJjnjfiie  I'roper- 
t-ica  of  Hi^hi-o-roniscJ  Aiv,  SWC-TR-56-35,  National  Bureau  of  Standards, 
Washington,  DC,  April  1957. 

11.  Doan,  L.  R.,  and  Nickel,  G.  H.,  A  :>ubvout  ine  fov  the  t  at  Ion  of  State  of 

Aiv ,  RTD  (WLR)  TM-63-2,  Air  Force  Weapons  Laboratory,  Kirtland  Air  Force 
Base,  New  Mexico,  May  1963. 


for  plots  of  histories  of  pressure,  pressure-impulse,  density,  velocity, 
dynamic  pressure,  and  dynamic  pressure-impulse. 

Further,  a  provision  for  energy  loss  due  to  radiation  is  provided  with  three 
different  techniques.  And  finally,  to  provide  a  sharply  defined  stable 
shock  front,  there  are  three  viscosity  options:  1)  no  artificial  viscosity 
used,  2)  addition  of  a  constant  viscosity  term,  and  3)  the  addition  of 
a  pseudoviscous  pressure  (Ref.  12). 

HULL  HYDROCODE 

HULL  is  a  2-D  hydrocode  written  for  either  Cartesian  or  cylindrical  geometries 
(Ref.  9).  It  is  an  outgrowth  of  the  SHELL-OIL  code  developed  by  Johnson  in 
the  early  1960s.  Matuska  and  Durrett  made  the  code  second-order  accurate  in 
time  and  space  in  the  Lagrangian  phase  and  developed  the  SAIL  software  manage¬ 
ment  system  in  1971  (Ref.  8),  thus  providing  accuracy  and  flexibility  not 
previously  possible.  A  3-D  Cartesian  version  is  under  development. 

The  HULL  code  was  originally  designed  to  model  fluid  behavior  in  a  multi¬ 
dimensional  Eulerian  continuum.  As  a  result  of  an  extensive  and  on-going 
development  program  this  code  has  evolved  into  a  sophisticated  and  versatile 
computational  tool.  Among  its  capabilities  are  the  simulation  of  high-explo¬ 
sive  detonations,  nuclear  weapons  effects  (including  airblast  precursors), 
diffusion  limitation  for  multimaterial  environments,  and  so  on.  HULL  is  used 
in  conjunction  with  the  SAIL  preprocessor  program  which  efficiently  tailors 
the  calculational  coding  to  the  problem  specifications. 

The  basic  HULL  code  solves  finite-difference  analogs  for  a  set  of  partial 
differential  equations  which  govern  the  behavior  of  a  compressible,  noncon¬ 
ducting,  inviscid  fluid.  The  local  state  variables  are  updated  in  two  phases. 

12.  Needham,  C.  E.,  I'cvclo; mcnt  of  an  Artificial  Viscosity  Function,  AFWL-TR- 
77-53,  Air  Force  Weapons  Laboratory,  Kirtland  Air  Force  Base,  New  Mexico, 
August  1977. 


First,  the  velocity  and  energy  state  of  the  fluid  are  advanced  in  time  by  a 
Lagrangian  calculation.  In  the  absence  of  mesh  gradients  this  method  is 
fully  second-order  accurate  in  the  Lagrangian  phase.  Then,  instead  of 
reconfiguring  the  calculational  mesh  to  complete  the  Lagrangian  phase,  final 
fluid  properties  are  redefined  in  an  Eulerian  reference  system  by  taking 
into  account  the  flux  of  mass,  momentum,  and  energy. 

Two  features  of  HULL  are  of  particular  importance  when  considering  REST  and/or 
DABS  simulation  experiments.  Detonation  phenomena  can  be  simulated  by 
routines  which  both  model  the  burning  of  an  explosive  and  limit  diffusion  when 
treating  several  material  species  in  the  Eulerian  framework.  An  algorithm 
controls  the  burning  process  so  that  the  detonation  propagates  sequentially 
through  the  mesh  at  the  local  sound  speed  and  in  the  direction  of  the  local 
pressure  gradient.  Equations  of  state  for  both  the  unburned  explosive  and 
detonation  products  are  used  to  update  the  relative  masses  and  the  internal 
energy  as  the  detonation  progresses.  The  HULL  code  diffusion  limiter  plays  a 
significant  role  in  such  a  multimaterial  environment.  It  tends  to  restrict 
the  mixing  of  specie s  in  the  direction  of  flow  and  thereby  preserves  contact 
discontinuities  which  are  characteristic  of  explosively  generated  blast 
waves.  Diffusion  is  constrained  by  an  algorithm  which  arbitrarily  adjusts  the 
flux  of  mass  by  giving  preference  to  those  species  already  downstream. 

Options  available  for  use  with  HULL  include  multimaterial  modeling,  high- 
explosive  burn,  immovable  islands,  stations,  and  special  input-boundary 
conditions.  The  islands  permit  calculation  of  shock  interactions  with 
structures.  The  boundary  conditions  permit  input  from  previous  SAP  or  HULL 
calculations.  The  stations  allow  description  of  all  the  hydrodynamic 
variables  as  a  function  of  time  at  any  designated  point  in  the  calculational 
space. 
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SECTION  III 

HEST  SIMULATOR  CALCULATIONS 


The  pressures  and  shock  oscillations  produced  during  the  explosion  of  a 
HEST  cavity  are  of  substantial  interest  to  the  community  of  research  in  simu¬ 
lation  development.  A  number  of  complex  hydrocode  calculations  have  been 
performed  to  further  understand  and  define  these  HEST  characteristics.  To 
achieve  a  crude  approximation  of  HEST  cavity  pressures,  a  1-D  similarity 
solution  was  generated  for  comparison  with  a  recently  performed  HEST  experi¬ 
ment.  The  HEST  Over  Rectangular  Slab  (HORS)  1-3  experiment  (Ref.  13) 
produced  good  pressure  data  and  was  selected  for  calculational  modeling  using 
a  1-D  similarity  approximation.  The  HORS  1-3  HEST  configuration,  depicted 
in  Figure  2,  contains  five  planes  of  explosive  materials  in  a  355.6-mm  explo¬ 
sion  cavity.  A  first  approximation  for  this  experiment  is  achieved  by  consi¬ 
dering  the  five  explosive  planes  in  the  experiment  to  be  consolidated  into  a 
single  plane  located  along  the  midplane  of  the  cavity  (Fig.  3).  The  charge 
density  of  2  x  10'r  g/mm3  used  in  the  experiment  is  equivalent  to  a  planar 
charge  density  of  4.297  x  108  ergs / mm2 .  The  foam-filled  cavity  used  in  HORS 
1-3  is  considered  to  be  filled  with  an  ideal  gas  (y  =  1.4).  Upon  initiation, 
shock  waves  propagate  from  either  side  of  the  explosive  plane  and  eventually 
impinge  on  the  cavity  boundaries.  Since  the  cavity  is  symmetrical  about  the 
explosive  plane,  only  a  half-space  need  be  calculated. 


ONE-DIMENSIONAL  SIMILARITY  S0LUTI0N--HEST  SINGLE  EXPLOSIVE  PLANE 

The  results  from  the  similarity  solution  are  shown  in  Figure  4.  To  carry 
the  solution  beyond  the  point  of  arrival  of  shock  at  the  boundary  is  not 
valid  because  the  similarity  solution  does  not  hold  for  reflected  shocks. 


13.  Gagnon,  L.  W.,  HEST  One >’  Rcatanjulae  Slab  (HORS.l—Phasc  7,  AFWL-TR-78-238, 
Air  Force  Weapons  Laboratory,  Kirtland  Air  Force  Base,  New  Mexico,  March 
1  979. 


Figure  2.  HORS  1-3  HEST  Experiment  Configuration 
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Figure  3.  One-Dimensional  Approximation  of  HEST  Explosion  Cavity 
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Figure  4.  Similarity  Solution--Planar  Blast  Wave 


The  solution  shows  the  shock  arrival  at  the  cavity  edge  (177.8  mm)  at  11.87 
with  a  peak  pressure  of  83.5  MPa  which  would  cause  reflected  shocks  in  the 
neighborhood  of  665  MPa  (Ref.  14).  The  cavity  equilibrium  pressure  near 
the  center  of  the  cavity  at  this  time  is  approximately  32  MPa.  The  data 
from  the  experiment  showed  the  cavity  equilibrium  pressure  (peak  simulation 
pressure)  to  be  approximately  34.65  MPa  and  the  peak  recorded  pressure  spikes 
to  be  in  excess  of  207  MPa.  There  is  some  question  as  to  whether  the  instru¬ 
mentation  system  was  able  to  respond  to  higher  frequency  transient  pressure 
spikes.  Hence  reflected  shocks  with  magnitudes  in  the  range  of  665  MPa 
may  have  occurred  but  were  not  recorded. 

A  second  1-D  model  of  the  HEST  cavity  contained  five  explosive  planes  as 
shown  in  Figure  2.  The  same  amount  of  explosive  energy  as  in  the  previous 
problem  (single  midheight  plane--  4.297  x  10p  ergs/mnr)  was  considered  to 
be  distributed  uniformly  among  the  five  planes  of  explosive  material,  i.e., 
8.594  x  107  ergs/mm'  at  each  plane.  The  explosion  from  the  plane  nearest 
the  cavity  boundary  was  calculated  out  to  the  time  of  arrival  at  the  cavity 
boundary  using  the  1-D  similarity  solution.  The  results,  plotted  in 
Figure  5,  show  the  shock  arriving  at  the  cavity  boundary  at  5.12  Hs  with 
a  peak  incident  pressure  of  49.728  MPa  which  is  substantially  less  than 
the  83.5  MPa  obtained  with  the  single  midheight  plane  of  explosives.  Hence 
it  appears  that  a  greater  dispersion  of  the  explosives  in  the  HEST  cavity 
leads  to  a  reduction  in  the  peak  incident  pressure  arriving  at  the  cavity 
boundary  (and  a  reduction  in  resulting  reflected  pressures)  along  with  an 
earlier  arrival  of  the  first  shock. 

The  1-D  similarity  solution  has  provided  a  first-order  approximation  for  the 
initial  blast  wave  produced  in  HEST  experiments  which  can  serve  as  a  basis 
for  comparison  for  the  hydrocode  calculations. 


Explosiv e  Shocks  in  Air ,  McMillan,  p.  57.,  1962 


14.  Kinney,  G.F., 
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Figure  5.  One-Dimensional  Similarity  Solution  (Planar  Blast  Wave-- 
Multiple  Explosive  Planes) 
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HULL  HYDROCODE  CALCULATION- -HEST  SINGLE  EXPLOSIVE  PLANE* 


A  HULL  calculation  (No.  905.0020)**  was  performed  for  comparison  with  the  1-D 
similarity  solution  discussed  previously  (Fig.  3).  A  6-  by  72-cell  Cartesian 
mesh  was  used  with  the  left  boundary  considered  to  be  the  midheight  of 
the  HEST  cavity  (Fig.  6).  The  right  boundary  was  set  to  be  reflective  in 
order  to  simulate  the  floor  of  the  HEST  explosion  cavity.  Top  and  bottom 
boundaries  were  also  set  to  be  reflective.  The  mesh  was  loaded  initially 
with  an  ideal  gas  having  a  gamma  of  1.4  and  a  density  of  1.225  x  10”s  g/mm3. 

A  quantity  of  energy  (5.4567  x  108  ergs)  was  introduced  into  each  of  the  cells 
along  the  left  boundary  at  time  zero  to  initiate  the  explosion  process.  This 
deposition  resulted  in  a  planar  energy  density  along  the  midplane  of  4.358  x 
108  ergs/mm2  (half  being  introduced  on  each  side  of  the  midplane).  The  value 
of  4.358  x  108  differs  slightly  from  the  4.297  x  108  value  used  in  the  1-D 
similarity  solution  because  of  an  adjustment  in  the  mesh  cell  size  which  was 
made  while  setting  up  the  problem. 

The  results  of  this  HULL  calculation  are  shown  in  Figures  7  through  9. 

Figure  7a  shows  the  calculated  shock  wave  at  2  ys  moving  from  the 
plane  of  symmetry  (left  boundary)  toward  the  right  boundary  (located  at 
177.8  mm).  Figure  7b  shows  the  shock  immediately  before  arrival  at  the 
reflecting  boundary.  The  time  is  12  ys  and  the  peak  overpressure  is  81.31  MPa, 
which  agrees  well  with  the  11.87  ys  and  83.87  MPa  obtained  from  the  similarity 
solution. 

Shock  reflections  from  rigid  surfaces  can  be  calculated  with  the  HULL  code. 
Figure  8a  shows  the  calculated  wave  profile  at  16  ys,  a  short  time  after 
reflecting  at  the  right  boundary.  It  is  shown  to  be  moving  back  toward  the 
plane  of  symmetry  and  has  a  peak  pressure  of  389  MPa.  Figure  8b  shows  the 
profile  of  the  reflected  wave  at  21  ys  as  it  continues  to  move  toward  the 
left  with  the  peak  pressure  now  down  to  114  MPa. 


*A  table  listing  all  of  the  hydrocode  calculation  data  reports  produced 
during  this  technical  effort  is  provided  in  Appendix  B. 

**Hydrocode  Calculation  905.0020,  Data  Report,  UNM/CERF  AST-35,  December  1978. 
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Pressure,  MPa 


Figures  9a  and  9b  are  history  plots  of  calculated  pressures  for  two  gage 
stations.  Station  38  is  at  the  right  reflecting  boundary  (177.8  mm)  and 
Station  18  is  midway  between  the  plane  of  symmetry  and  the  right  boundary.  The 
peak  reflected  pressure  at  the  right  boundary  is  589  MPa,  ocurring  at  14.57  ps. 
The  readings  are  in  reasonable  agreement  with  the  estimate  of  665  MPa  based 
on  the  similarity  solution.  The  arrival  and  reflection  of  the  shock  from 
the  right  boundary  a  second  time  can  also  be  seen  in  Figure  9a.  The  peak 
reflected  pressure  is  substantially  lower  the  second  time  (140  MPa)  due  to 
the  second  shock  traveling  through  previously  shocked  ideal  gas  and  to 
pseuaoviscous  effects  innerent  in  the  HULL  code.  Figure  9b  shows  four  distinct 
shock  arrivals  which  are  identified  as  1)  the  first  arrival  of  the  incident 
wave  traveling  from  the  plane  of  symmetry  toward  the  right  boundary, 

2)  the  arrival  of  the  reflected  wave  from  the  right  boundary,  3)  the  arrival 
of  the  reflected  wave  from  the  plane  of  symmetry,  and  4)  the  arrival  of 
the  second  reflection  from  the  right  boundary.  The  small  oscillations  in  the 
waveform  following  the  third  and  fourth  peaks  are  caused  by  instabilities 
in  the  calculation. 

From  these  results  it  can  be  seen  that  in  a  HEST  experiment  where  a  single 
plane  of  explosive  is  placed  at  the  midheight  of  the  explosion  cavity,  a 
large  pressure  spike,  many  times  the  magnitude  of  the  cavity  equilibrium 
pressure,  is  to  be  expected  at  the  floor  of  the  cavity.  Numerous  shock 
oscillations  can  be  expected  throughout  the  cavity  during  the  first  few 
milliseconds  of  the  experiment. 

Calculation  905.0020  was  repeated  in  calculation  905.0050*  except  that  the 
equation  of  state  [P  =  f(p,  E)]  described  the  explosion  cavity  gases  as 
real  air  instead  of  an  ideal  gas.  The  resulting  pressure  history  plots  for 
Stations  38  and  18  (which  are  the  same  stations  shown  in  Figure  9)  are 
shown  in  Figure  10.  Although  the  wave  shapes  are  similar  a  significant 
reduction  in  peak  pressures  (approximately  one-half)  and  a  delay  in  arrival 
times  can  be  observed  when  compared  with  Figure  9.  This  would  tend  to 

♦Hydrocode  Calculation  905.0050,  Data  Report,  UNM/CERF  AST-33,  December  1978. 


discredit  the  previously  mentioned  reasonable  agreement  between  the  1-D 
similarity  solution  and  the  HORS  1-3  experiment.  Such  may  not  be  the  case. 
Recalling  that  the  explosion  cavity  in  the  HORS  1-3  experiment  was  filled 
with  a  polystyrene  plastic  foam  material,  it  is  quite  possible  that  the 
departure  from  real  air  behavior  caused  by  the  plastic  foam  would  cause 
similar  results  to  those  caused  by  the  use  of  an  ideal  gas  EOS  (e.g.,  higher 
peak  pressures  and  more  rapid  shock  propagation).  This  issue  is  not  fully 
resolved  and  more  study  is  required  to  completely  answer  the  question. 

ONE-DIMENSIONAL  HYDROCODE  CALCULATIONS --HEST  SINGLE  EXPLOSIVE  PLANE 

The  objective  of  this  calculational  study  was  to  evaluate  various  other  methods 
for  computer  simulation  of  the  test  environment  existing  in  a  HEST  explosion 
cavity.  Calculations  were  performed  using  the  SAP  and  APOD  hydrocodes. 

The  configuration  selected  as  a  baseline  for  comparison  was  an  air-filled 
cavity  as  shown  in  Figure  3.  The  plane  of  symmetry  represents  the  cavity 
midheight  as  well  as  the  location  of  the  explosive  plane.  The  specifications 
for  the  HORS  1-3  experiment  (355.6  mm  cavity  height  and  4.297  x  10B  ergs/mm2 
planar  explosive  charge  density)  were  used  to  model  the  calculations  within  a 
Cartesian  mesh.  Detonation  of  the  explosive  plane  was  modeled  to  occur 
instantaneously  along  this  midplane. 

Five  calculations  (905.0021,  905.0022,  905.0031,  905.0041,  and  905.0051)*  were 
made  simulating  the  one-dimensional  flow  between  the  plane  of  symmetry  and 
the  reflecting  soil  boundary.  The  APOD  hydrocode  was  used  for  905.0022 
and  the  SAP  hydrocode  was  used  for  the  rest.  The  plane  of  symmetry  was 
also  treated  as  a  reflecting  boundary.  Four  of  the  calculations  modeled 
the  material  within  the  main  cavity  as  an  ideal  gas  with  a  constant  value 
for  gamma  (1.4).  The  fifth  calculation  modeled  this  region  with  a  real 


*Hydrocode  Calculation  905.0021,  Data  Report,  UNM/CERF  AST-34,  December  1978 
Hydrocode  Calculation  905.0022,  Data  Report,  UNM/CERF  AST-22,  December  1978. 
Hydrocode  Calculation  905.0031,  Data  Report,  UNM/CERF  AST-21,  December  1978. 
Hydrocode  Calculation  905.0041,  Data  Report,  UNM/CERF  AST-23,  December  1978. 
Hydrocode  Calculation  905.0051,  Data  Report,  UNM/CERF  AST-24,  December  1978. 
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Figure  11.  One-Dimensional  Model  of  HEST  Cavity 

air  EOS.  The  mesh  was  initialized  with  a  mass  density  of  1.225  x  10*sg/mm3 
and  an  energy  density  of  2.067484  x  1 09  ergs/g  which  yielded  an  ambient 
pressure  of  0.102366  MPa  within  the  reaction  cavity.  Energy,  corresponding 
to  one-half  the  charge  density,  was  deposited  in  the  cell  next  to  the  plane 
of  symmetry  as  shown  in  Figure  11. 

Data  history  stations  were  placed  as  shown  in  Table  1  for  the  three  hydrocodes 
used  (the  HULL  calculation  905.0020  has  been  included).  Station  radius  is 
the  distance  of  the  station  from  the  plane  of  symmetry.  The  reflecting 
boundary  was  located  at  177.8  mm  (between  SAP  Stations  38  and  39).  Stations 
beyond  177.8  mm  were  used  in  SAP  calculations  905.0031  and  905.0041  which 
treated  the  reflecting  boundary  as  the  surface  of  a  very  dense  gas  rather 
than  a  rigid  wall , 

Calculations  905.0020,  905.0021,  and  905.0022  were  run  to  formulate  a  direct 
comparison  of  the  results  of  the  HULL,  SAP,  and  APOD  hydrocodes  respectively. 
Input  parameters  were  specified  to  make  the  initial  conditions  in  the  three 
calculations  agree  as  closely  as  possible.  APOD  is  a  purely  Lagrangian  code 
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while  both  HULL  and  SAP  contain  both  Eulerian  and  Lagrangian  features. 

Figures  12,  13,  and  14  show  comparable  pressure  history  plots  at  two  stations 
for  each  of  the  three  calculations. 


Table  2  presents  a  comparison  of  the  peak  values  at  these  stations.  The 
lower  curve  of  Figure  12  represents  the  HULL  solution  for  the  station 
immediately  in  front  of  the  reflecting  boundary  (Station  38).  The  peak 
reflected  pressure  at  the  right  boundary  is  589.84  MPa  occurring  at  14.577  ys . 
The  reflection  of  the  shock  from  the  right  boundary  a  second  time  can  also 
be  noted.  This  second  peak  is  substantially  lower  (140  MPa)  due  to  the 
change  in  properties  of  the  ideal  gas  after  the  passage  of  the  first  shock  and 
to  psuedoviscous  effects  within  the  HULL  hydrocode.  The  upper  curve  of 
Figure  12  is  for  Station  30  which  is  141.48  mm  from  the  plane  of  symmetry 
and  shows  four  distinct  shock  arrivals  which  are  identified  from  left  to 
right  as  1)  the  initial  arrival  of  the  incident  wave  traveling  from  left 
to  right,  2)  the  arrival  of  the  reflected  wave  from  the  right  boundary, 

3)  the  arrival  of  a  wave  reflected  from  the  plane  of  symmetry,  and  4)  the 
arrival  of  the  second  reflection  from  the  right  boundary.  The  HULL  and  SAP 
curves  (Figs.  12  and  13)  are  almost  identical  in  shape  and  arrival  time. 

The  APOD  curves  (Fig.  14)  have  the  same  general  form  but  the  peak  values 
are  somewhat  lower  and  the  shock  waveforms  are  definitely  smeared,  taking 
on  an  overdamped  appearance  rather  than  the  sharp  definition  of  the  other 
two  solutions. 

In  calculations  905.0031  and  905.0041  the  reflecting  wall  of  the  HEST  cavity 
was  simulated  with  a  nonrigid,  partially  reflecting  heavy  gas  (the  rigid  wall 
of  905.0021  was  replaced  with  additional  cells  of  ideal  gas  which  had  a  mass 
density  of  1.602  x  10"?  g/mm3).  A  gamma  coefficient  of  1.05  and  an  energy 
density  of  1.279  x  107  ergs/g  were  selected  for  905.0031,  while  905.0041  used 
the  same  gamma  (1.4)  used  in  the  main  cavity  and  an  energy  density  of 
1.581  x  106  ergs/g.  In  both  calculations  the  energy  density  was  computed 
as  a  function  of  the  specified  mass  density  and  gamma  so  that  the  ambient 
pressure  would  match  that  of  the  main  cavity.  Figures  15  and  16  show  the 
pressure  histories  for  these  two  calculations  at  the  same  two  stations  as 
shown  previously  (Figs.  12,  13,  and  14).  Note  in  Table  2  and  Figure  15  that 
calculation  905.0041  produces  the  first  reflections  with  reduced  peaks 
and  slightly  delayed  arrival  times,  but  that  all  of  the  secondary  reflections 
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.  SAP  Calculation  905.0021 --Pressure  versus  Time 
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are  quite  different  with  much  later  arrival  times.  However,  calculation 
905.0031  is  an  entirely  different  story;  here  the  first  reflection  is 
extremely  weak.  Note  also  in  Table  2  and  Figure  15  at  141.48  mm  the 
reflected  peak  is  lower  than  the  incident  peak,  contrary  to  all  of  the  other 
calculations.  This  method  of  modeling  a  nonrigid  reflecting  boundary  does 
not  appear  to  t  satisfactory. 

Calculation  905.0051  is  a  modification  of  905.0021  with  the  ideal  gas  re¬ 
placed  with  real  air  employing  a  variable  gamma.  The  reflecting  surface 
represents  rigid,  fully  reflecting  wall.  Even  with  the  same  initial  mass  and 
energy  deposition,  the  equation  of  state  employed  in  905.0051  produces  an 
initial  peak  pressure  of  1807.1  MPa  which  is  approximately  one-half  the 
comparable  peak  of  3473.6  MPa  for  calculation  905.0021.  This  trend  is 
evident  at  the  other  stations  (Table  2  and  Fig.  17).  The  arrival  times  in 
the  real  air  case  were  considerably  later  than  those  for  the  ideal  gas. 

These  results  are  consistent  with  the  earlier  findings  in  calculations 
905.0020  and  905.0050. 

Computer  costs  for  performing  SAP,  APOD,  and  HULL  calculations  were  obtained 
during  this  study.  The  SAP  and  APOD  calculation  costs  were  approximately  the 
same,  $10.50  for  APOD  and  $8.00  for  SAP.  The  HULL  calculation  cost  $95.00. 

All  calculations  were  run  to  80  ;.s. 

The  APOD  hydrocode  was  dropped  from  any  further  utilization  during  this  techni¬ 
cal  effort.  The  versatility  of  SAP,  having  both  a  Lagrangian  and  Eulerian 
capability,  as  well  as  an  excellent  supporting  executive  program  and  a  strong 
plot  program  at  the  same  cost  level  as  APOD,  makes  SAP  far  more  attractive 
for  one-dimensional  problems.  Additionally  SAP  provided  results  almost 
identical  to  HULL  results.  The  one  advantage  held  by  APOD  is  that  of  having 
a  1  1/2-D  capability  which  will  be  overcome  with  the  proposed  SAP  1  1/2-D 
modification.  SAP  will  then  be  used  for  one-dimensional  problems  and  certain 
1  1/2-D  problems  while  HULL  will  be  used  in  2-D  applications. 

The  work  described  above  has  clearly  shown  the  capability  of  the  SAP  and  HULL 
codes  to  qualitatively  model  a  number  of  1-D  shock  reflection  and  multiple 


Pressure,  MPa  Pressure,  MPa 


shock  problems  encountered  in  HEST  simulation  experiments.  Additional 
work  is  now  required  to  validate  these  hydrocodes  quantitatively  using 
high  quality  experimental  data. 

HULL  SIMULATIONS  OF  DETONATING  CORD 
BLAST  EFFECTS  IN  A  HEST  CAVITY 

The  explosive  charge  used  with  the  HEST  simulator  has  generally  been  fabri¬ 
cated  using  individual  strands  of  detonating  cord.  A  matter  of  concern  to 
HEST  designers  lias  been  the  severity  of  the  blast  environment  to  which 
pressure  gages  are  exposed.  The  validity  of  data  obtained  from  gages  in  the 
vicinity  of  detonating  cords  is  frequently  questioned.  To  better  understand 
and  help  resolve  this  issue  a  series  of  four  HULL  calculations  was  performed. 
Individual  strands  of  detonating  cord  were  coarsely  modeled  in  layered  planes 
within  a  typical  HEST  explosion  cavity.  Of  particular  interest  were  the 
reflected  peak  pressures  experienced  at  the  cavity  floor  where  gages  are  often 
installed. 

The  HEST  explosion  cavity  was  assumed  to  be  initially  filled  with  an  ideal 
gas  (y  =  1.4)  having  properties  approximating  a  standard  sea  level  atmosphere. 
All  boundaries  of  the  calculational  mesh  were  specified  as  reflective  to 
represent  the  cavity  walls  as  fixed  surfaces.  For  the  purposes  of  these 
calculations  the  cavity  was  assumed  to  contain  five  equally  spaced  horizontal 
explosive  planes  composed  of  equally  spaced  strands  of  detonating  cord 
(Fig.  18).  The  lateral  distance  between  strands  was  taken  to  be  equal  to  the 
distance  between  planes. 

Variations  of  the  basic  detonating  cord  pattern  were  simulated  in  a  2-D 
Cartesian  geometry  where  the  plane  of  symmetry  represented  the  cavity 
midheight.  Individual  strands  were  represented  with  rectangular  blocks  of 
20  cells  so  that  the  cross-sectional  area  corresponded  to  a  cord  diameter  of 
13  mm. 

Detonations  were  assumed  to  be  instantaneous  as  well  as  simultaneous  and  were 
simulated  by  depositing  energy  into  the  designated  cells  between  the  second 
and  third  calculational  cycles.  In  all  cases  this  deposition  was  apportionea 


equally  among  the  explosive  planes  to  achieve  an  overall  horizontal  planar 
density  similar  to  that  of  the  HORS  1-3  experiment  (4.297  x  108  ergs/mm2). 


This  calculational  study  must  be  considered  qualitative  since  several  physical 
phenomena  were  not  modeled  in  the  calculations.  Not  considered  were  the  mass 
of  the  explosive  material,  the  detonation  process,  and  the  products  of  combus¬ 
tion.  Consequently  detonations  are  simulated  by  the  rapid  expansion  of  a  very 
hot  ideal  gas.  A  baseline  calculation  (905.0030)  simulated  the  blast  environ¬ 
ment  produced  by  the  basic  detonating  cord  pattern.  In  succeeding  calculations 
(905.0040,  905.0060,  and  905.0090)  this  basic  pattern  was  variously  modified 
(Fig.  19)  in  attempts  to  locally  attenuate  the  peak  reflected  pressure  at 
the  right-hand  boundary  which  represents  the  cavity  floor  where  pressure 
measuring  gages  would  be  located. 


Calculation  905.0030* 

Fluid  properties  were  monitored  at  two  locations  on  the  right-hand 
reflective  boundary.  At  locations  immediately  opposite  a  detonating  cord 
(e.g..  Station  10)  the  calculated  peak  reflected  pressure  of  158  MPa  was 
caused  by  the  regular  reflection  of  the  shock  front  impinging  at  near  normal 
incidence.  The  maximum  reflected  peak  (227  MPa)  was  observed  at  the  floor 
location  approximately  midway  between  two  adjacent  detonating  cords  (Station  7). 
In  this  instance  it  is  believed  that  the  incident  wave  system  (Fig.  20)  is  a 
Mach  configuration  generated  by  the  simultaneous  explosion  of  two  equal  charges 
(Ref.  15).  The  plane  of  symmetry  midway  between  these  detonations  is 
essentially  a  reflective  surface.  Since  at  early  times  the  shocks  are  very 
strong,  nonregular  (i.e.,  Mach)  reflections  occur  and  give  rise  to  a  Mach 
shock  front  which  is  perpendicular  to  this  plane  and  which  is  generated  by 
the  interaction  of  the  two  cyl indrical ly  expanding  shock  waves  from  the 
individual  detonating  cords.  The  reflection  at  Station  7  is  most  likely  the 
result  of  the  normal  incidence  of  this  Mach  shock  front. 


15.  Courant  and  Friedrichs,  ,  rr.wA*  F  .mJ  -v.-j ,  pp.  334-5 

Interscience  Publishers,  New  York,  1948. 

*Hydrocode  Calculations  No.  905.0030  and  905.0070,  Data  Report,  UNM/CERF 
AST-43,  December  1978. 
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Figure  19.  Detonating  Cord  Patterns 
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Figure  20.  Incident  Wave  System— Calculation  905.0030 
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Calculation  905.0040* 


In  this  model  the  basic  pattern  is  interrupted  by  omitting  the  detonating  cord 
directly  above  Station  10.  Otherwise  the  model  specifications  are  identical 
to  those  for  calculation  905.0030.  An  illustration  of  the  wave  system  produced 
by  this  pattern  is  depicted  in  Figure  21.  As  this  wave  system  evolves  it  is 
symmetrical  with  respect  to  the  X-Z  plane  containing  Station  10.  The  reflected 
peak  pressure  of  91  MPa  at  Station  7  was  the  result  of  the  oblique  incidence 
of  the  shock  front  generated  by  the  detonation  of  the  cord  at  Ei.  The  634  MPa 
reflected  peak  pressure  observed  at  Station  10  involves  the  staging  of  several 
Mach  interactions  and  then  the  nearly  simultaneous  convergence  of  the  result¬ 
ing  ,'ct  and  the  reflections  of  the  two  oblique  shocks  at  the  cavity 
floor  (Fig.  22). 

The  first  shock  interactions  occur  midway  between  cords  Ej  and  E2,  and  E2  and 
E 3 .  The  Mi  shock  fronts  produced  here  converge  at  Fi  and  act  to  reinforce 
the  M2  shock  in  the  direction  of  F2.  The  convergence  of  the  M2  and  Mj  fronts 
at  F2  eventually  leads  to  the  formation  of  a  Mach  jet  since  the  flow  is 
locally  constrained  by  the  primary  shock  fronts  propagating  inward  from 
E3  and  E3.  The  evolution  of  this  wave  system  and  eventual  convergence  of 
shocks  at  F3  is  clearly  discernible  in  the  sequential  display  of  the  contoured 
pressure  fields  (Fig.  23). 

Calculation  905.0060** 

This  calculation  simulated  a  proposed  technique  for  suppressing  the  shock 
focusing  experienced  in  calculation  905.0040.  The  region  of  energy  deposition 
for  the  cord  immediately  above  Station  10  was  expanded  so  that  the  energy 
previously  deposited  in  a  few  cells  was  now  uniformly  distributed  among  506 
cells  (called  a  buffer  zone).  Otherwise  the  model  specifications  remained 
unchanged.  Since  the  initial  pressure  in  this  region  is  approximately  60  MPa, 
it  acts  as  a  local  buffer  zone  and  the  wave  system  seen  in  the  905.0040 


*Hydrocode  Calculaton  No.  905.0040,  Data  Report,  UNM/CERF  AST-39,  December 
1978. 

**Hydrocode  Calculation  No.  905.0060,  Data  Report,  UNM/CERF  AST-40,  December 
1978. 
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Figure  21.  Wave  System  at  4  ys--Calculation  905.C040 
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Calculation  905.0040--HULL  Simulation  of  HORS  1-3  HEST  Detonating 
Cord  Array  (3  of  10) 
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Figure  24.  Shock  Wave  Interactions--Calculation  905.0060 
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solution  did  not  fully  develop.  Instead  a  secondary  wave  system  developed 
within  this  expanding  zone  (Fig.  24)  and  the  reflected  peak  pressure  at  Station 
10  was  reduced  to  195  MPa.  The  initial  front  interaction  with  this  configur¬ 
ation  would  be  the  head-on  collisions  of  the  nearly  planar  Si  fronts  with  the 
much  stronger  S2 ,  S2,  and  S3  shocks.  These  shock  waves  penetrate  each  other 
and  in  the  process  both  are  attenuated.  At  approximately  4  us  the  Si  shock 
wave  impinges  on  the  cavity  floor  at  Station  10.  Meanwhile  the  weakened  S2 
and  S )  shocks  propagate  into  the  buffer  zone  and  interact  at  the  X-Z  plane 
containing  Station  10.  This  plane  functions  as  a  reflecting  surface,  and  the 
regular  reflections  of  S2  and  S2  eventually  act  to  retard  the  Si  reflection 
advancing  upward  from  Station  10.  The  constraint  of  this  converging  wave 
system  may  be  responsible  for  the  pressurization  at  this  location.  Finally  the 
reflected  peak  pressure  observed  at  Station  7  (221  MPa)  is  believed  to  be 
caused  by  the  convergence  of  the  regular  reflections  of  the  primary  Si  and  S' 
waves.  Although  these  results  do  not  compare  too  favorably  with  those  produced 
by  the  original  basic  pattern  (905.0030),  they  do  demonstrate  the  feasibility 
of  the  protective  buffering  technique. 

Calculation  9C5.0090* 

This  calculation  is  actually  a  sequel  to  calculation  905.0060.  To  evaluate 
its  attenuating  capability  the  buffer  zone  was  expanded  to  include  2001  cells 
among  which  energy  equivalent  to  the  deposition  for  3  detonating  cords  was 
evenly  distributed.  As  a  result  the  initial  pressure  in  this  zone  was  reduced 
to  46.5  MPa.  Since,  in  this  case,  the  lower  boundary  of  the  zone  is  approx¬ 
imately  10  mm  from  the  cavity  floor,  the  initial  reflected  peak  pressure  at 
both  Stations  10  and  7  (80  MPa  at  3.5  s)  is  simply  due  to  the  early  time 
expansion  of  this  buffer  zone.  However,  late-time  reflections  are  substan¬ 
tially  greater.  Peaks  of  99  MPa  at  12.5  ys  and  108  MPa  at  13.5  ys  are 
experienced  at  Stations  7  and  10  respectively. 

The  evolution  of  the  pressure  field  leading  up  to  these  secondary  reflections 
is  not  clear,  but  some  general  observations  can  be  made.  The  complex  wave 
system  spawned  by  the  shock  interactions  within  the  cavity  interior  undergoes 

*  Hydrocode  Calculation  No.  905.0090,  Data  Report,  UNM/CERFAST-41 ,  December 
1  978. 
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a  two-stage  attenuation.  At  approximately  1.5  ps  the  upper  plane  wave  from 
the  buffer  zone  collides  with  the  convex  shock  fronts  produced  by  detonating 
cords  immediately  above.  As  these  weakened  and  retarded  fronts  descend 
through  the  buffer  zone  the  center  portion  of  the  wave  system  eventually 
develops  into  a  shallow  concave  front  followed  by  a  local  high-pressure 
wedge.  Meanwhile  the  reflected  planar  shock  front  has  been  propagating  upward 
from  the  cavity  floor  (the  pressure  histories  for  Stations  6  and  9*  would 
indicate  an  arrival  at  approximately  7  ps).  As  these  opposing  fronts  meet, 
the  shocks  approaching  the  cavity  floor  are  once  again  slowed  and  diminished. 
The  second  collision  apparently  occurs  at  9  ps  and  may  be  correlated  with  the 
high-pressure  ridge  which  is  observed  in  the  contoured  pressure  field  (Fig. 
25).  Once  these  attenuated  shocks  impinge  on  the  floor  the  subsequent 
interaction  of  the  reflections  is  quite  similar  to  that  observed  during 
calculation  905.0060. 

A  final  consideration  in  this  study  was  to  evaluate  the  case  where  the 
buffer  zone  represents  the  entire  cavity  (representing  cavity  equilibrium 
pressure).  Assuming  the  cavity  is  perfectly  insulated  and  maintains  its 
original  size  and  shape,  a  solution  of  the  EOS  would  yield  the  cquilibvi.c 
pressure  (i.e.,  the  uniform  pressure  within  the  cavity  after  an  infinite 
time).**  Given  that  for  an  ideal  polytropic  gas 

P  =  kpe  (y  -  1 )  (1 ) 

where 

P  =  pressure  (MPa) 

.  =  mass  density  =  1.225  x  10”5  g/mrrr 
=  specific  internal  energy  (ergs/g) 

,  =  adiabatic  exponent  =  1 .40457725 
k  =  conversion  factor  =  10~4 

it  is  first  necessary  to  solve  for  the  cavity  uniform  energy  density.  Given 
the  cavity  horizontal  planar  energy  density  (4.297  x  10p  ergs/mm'),  it  then 
follows  that 

♦Located  30  mm  above  the  cavity  floor. 

**The  equilibrium  pressure  correspords  most  closely  to  the  peak  simulation 
pressure  (PP  )  defined  in  Figure  1. 
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Figure  25.  Calculation  905.0090--HULL  Simulation  of  HORS  1-3  HEST  Detonating 
Cord  Array 
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4.297  x  1 0 8 / ( 1 4  x  25.4) 

L  ~  - - - * - 

1.225  x  10'6 

and  therefore  that 

P  =  [4.297  x  1 08/ ( 1 4  x  25.4)]  (y  -  DOO-1*) 

=48.9  MPa 

HULL  calculation  905.0080  was  run  with  the  explosive  energy  uniformly  intro¬ 
duced  in  each  of  the  cells  in  the  calculational  mesh  to  obtain  a  computed 
cavity  equilibrium  pressure.  The  calculation  was  run  twice,  once  with  an 
ideal  gas  and  once  with  a  real  air  EOS,  resulting  in  equilibrium  pressures  of 
48.99  MPa  and  20.48  MPa  respectively.  Since  the  calculation  was  run  for  only 
a  few  cycles,  no  data  reports  were  produced. 

A  comparison  of  the  calculation  results  and  the  equilibrium  pressure  is  shown 
in  Figures  26  through  28.  Calculation  905.0020  was  included  in  Figure  26 
to  provide  a  comparison  with  the  case  where  all  the  explosive  energy  was 
deposited  in  the  column  of  cells  adjacent  to  the  plane  of  symmetry  at  cavity 
midheight. 

These  calculations  clearly  demonstrate  that  complex  shock  interactions  which 
take  place  within  an  explosively  driven  test  cavity  can  result  in  severe  local 
pressure  concentrations.  This  is  particularly  true  for  the  typical  HEST  field 
experiment  where  the  desired  pressure  environment  is  generated  by  techniques 
using  concentrated  explosive  charges.  It  is  also  evident  that  it  may  be 
possible  to  distribute  the  explosive  energy  more  evenly  to  locally  attenuate 
the  pressure  field  at  gage  locations. 

However,  when  one  considers  applying  the  protective  buffering  technique  used 
in  calculations  905.0060  and  905.0090,  some  practical  constraints  must  be 
borne  in  mind.  An  absolutely  uniform  energy  distribution  throughout  such  a 
buffer  zone  cannot  be  achieved  experimentally.  There  are,  for  instance, 
obvious  limitations  if  one  elects  to  refine  the  detonating  cord  pattern  in 
order  to  produce  the  buffer.  Even  with  the  use  of  the  smallest  cord 
available,  reflected  peak  pressures  could  be  appreciable.  Other  seemingly 
more  appropriate  techniques  such  as  the  use  of  an  explosive  gaseous  mixture 
have  proven  unsatisfactory  because  of  safety  and  engineering  considerations. 
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(2) 

(3) 


Fiqure  26.  Peak  Reflected  Pressure  at  Cavity  Floor 


62 


Pressure,  MPa  Pressure,  MPa 


Equilibrium  Pressure 


Figure  27.  Station  7  Pressure 


Pressure,  MPa  Pressure,  MPa 
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In  summary  the  results  of  these  calculations  are  not  quantitatively  con- 
conlusive.  For  example,  the  aforementioned  inadequacies  in  the  simulation  of 
detonation  phenomena  may  have  had  a  marked  effect  on  these  solutions. 

To  relate  such  numerical  simulations  more  meaningfully  to  the  HEST  environment 
the  calculational  model  must  be  refined  in  several  respects. 

1.  The  material  model  for  the  detonating  cord  explosive  must  realisti¬ 
cally  represent: 

a.  the  mass  of  the  explosive  and  its  surrounding  plastic  case, 

b.  the  detonation  process,  and 

c.  the  products  of  detonation. 

2.  The  high  temperatures  encountered  in  blast  environments  require  that 
air  be  treated  as  a  real  gas  with  variable  specific  heats. 

3.  Protective  buffer  zones  should  be  simulated  as  the  result  of 
multiple  detonations  rather  than  as  a  uniform  pressure  field. 

Some  qualitative  conclusions  can  be  made.  The  practice  of  removing  strands 
of  detonating  cord  from  the  region  immediately  above  the  pressure-measuring 
location  should  be  discouraged.  Instead  the  gages  should  be  placed  directly 
in  line  beneath  the  detonating  cord  strand,  or  even  better,  placed  in  the  space 
midway  between  the  strands  and  the  line  between  the  strand  (midway  between 
Station  10  and  Station  7  in  the  preceding  calculations).  The  greatest 
disbursement  of  explosive  materials  in  the  explosion  cavity  should  be  used, 
consistent  with  cost  and  practical  construction  considerations.  The  range 
of  the  pressure-measuring  gage  should  be  specified,  based  on  calculated 
estimates  of  the  incident  and  reflected  peak  pressures  at  that  boundary  (much 
higher  than  the  anticipated  peak  simulation  pressure  for  the  experiment). 

HULL  SIMULATION  OF  HEST  SINGLE  EXPLOSIVE 
PLANE  WITH  SWEEPING  WAVE 

Calculations  905.0020,  905.0021,  and  905.0051  simulated  the  HORS  1-3  HEST  air 
cavity  with  the  simple  technique  of  a  one-time  energy  dump  into  a  1-D  mesh. 
Another  approach  to  modeling  this  experiment  was  attempted  utilizing  a  timed 
sequence  of  energy  dumps  along  the  plane  of  symmetry  of  a  2-D  mesh,  thus 
producing  a  sweeping  explosion  wave.  This  sweeping  wave  of  energy  dumps  more 
realistically  represents  the  exploding  edge  of  a  woven  layer  of  detonating 
cord  used  in  more  conventional  HEST  experiments. 
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In  Figure  29  the  Y-axis  represents  the  plane  of  symmetry  (centerline)  of  a 
HEST  explosion  cavity.  In  the  calculations  this  is  modeled  as  a  reflecting 
boundary  due  to  the  symmetry  of  the  cavity.  The  right  boundary  of  Figure  29 
is  also  a  reflecting  boundary  representing  the  soil  overburden.  Energy  is 
introduced  into  the  mesh  in  the  cells  adjacent  to  the  Y-axis,  starting  at  the 
origin  and  proceeding  in  the  positive  Y-direction  in  proportion  to  the 
computed  distance  traveled  at  a  specific  problem  time.  A  constant  detonation 
velocity  is  used.  The  position  of  the  sweeping  detonation  point  at  some 
arbitrary  time  is  represented  in  Figure  29. 


Reflecting  [ 
Boundary 
(Plane  of 
Symmetry) 

Sweeping  j 

Detonation  ^1 
Point  I 


177.8  mm 


Constant  Detonation 
Velocity 


Reflecting 

Boundary 

(Soil) 


Ori gi n 


■Midheight  of  Explosion  Cavity 


Figure  29.  HEST  Explosion  Cavity  with  Sweeping  Detonation 

Three  HULL  calculations  were  accomplished  using  the  sweeping  wave  simulation. 
In  each  of  these  calculations  a  cell  width  and  height  of  10.459  mm  was  used. 
The  right  reflecting  edge  was  located  at  177.8  mm.  The  mesh  was  considered 

to  be  air  having  a  mass  density  of  1.225  x  10"r'  g/mm?  and  an  energy  density 

of  2.067  x  109  ergs/g.  These  conditions  yielded  an  initial  cell  mass 

of  1.34  x  10~3  g,  an  initial  cell  energy  of  2.77  x  109  ergs,  and  an  ambient 

cavity  pressure  of  0.1025  MPa.  The  constant  velocity  of  detonation  was 
6.4008  x  109  mm/s.  A  real  air  EOS  was  used  with  a  variable  gamma  coefficient. 
Data  history  stations  were  placed  as  shown  in  Table  3  for  all  three  calcula¬ 
tions  . 


In  calculations  905.1080  and  906.1010*  the  energy  added  to  the  mesh  was 
5.459  x  109  ergs  per  cell  which  amounted  to  a  planar  charge  density  of  9.98  x 
107  ergs/mm2.  This  value  is  substantially  lower  than  the  desired  planar  charge 
density  of  4.297  x  10s  ergs/mm2  due  to  an  oversight.  In  both  calculations  the 
top  (sweep  exit)  boundary  was  transmissive. 


TABLE  3.  STATION  LOCATIONS 


Station 

No . 

Co-ordi nates 

X  Y 

mm 

Station 

No. 

Co-ordinates 

X  Y 

mm 

1 

170 

30 

29 

20 

180 

17 

170 

60 

30 

30 

180 

2 

170 

90 

31 

40 

180 

18 

20 

120 

32 

60 

180 

19 

30 

120 

33 

80 

18C 

20 

40 

120 

34 

90 

180 

21 

60 

120 

35 

110 

180 

22 

80 

120 

36 

130 

180 

23 

90 

120 

37 

150 

180 

24 

no 

120 

38 

170 

180 

25 

130 

120 

39 

170 

180 

26 

150 

120 

4 

170 

210 

27 

160 

120 

50 

170 

250 

28 

170 

120 

5 

170 

280 

3 

170 

150 

61 

170 

310 

Calculation  905.1080  utilized  a  transmissive  lower  (sweep  entry)  boundary 
while  calculation  906.1010  had  the  same  parameters  except  for  a  reflective 
lower  boundary.  Figure  30  shows  a  pressure  history  at  a  station  next  to  the 
right  reflecting  boundary  while  Figure  31  shows  a  plot  of  equal  pressure 
contours  and  a  plot  of  velocity  vectors  corresponding  approximately  to  the 
time  of  maximum  reflected  pressure  in  Figure  30.  It  appeared  that  the 
transmissive  lower  boundary  may  be  allowing  mass  and/or  energy  to  escape  from 
the  problem  mesh.  The  bottom  row  of  velocity  vectors  in  Figure  31  shows 
several  vectors  directed  outside  the  mesh  which  leads  to  this  interpretal 
Calculaton  906.1010  was  then  accomplished  with  the  lower  boundary  changed  to 


*Hydrocode  Calculation  No.  905.1080,  Data  Report,  UNM/CERF  AST-25,  December 
1  978. 

Hydrocode  Calculation  No.  906.1010.  Data  Report,  UNM/CERF  AST-29,  December 
1  978. 
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be  reflective  rather  than  transmitive.  The  reflective  lower  boundary  had  a 
very  significant  effect  on  the  total  results.  The  pressure  at  the  stations 
next  to  the  right  boundary  had  only  begun  to  show  a  pressure  increase  for  cal¬ 
culation  906.1010  when  the  problem  was  stopped  at  60  ps.  Note  that  in  Figure 
32  there  is  a  row  of  vectors  originating  from  a  lower  position  than  the  lowest 
row  of  Figure  31.  This  observation  leads  to  the  implication  that  the  trans¬ 
missive  boundary  of  calculation  905.1080  earlier  thought  to  be  allowing  mass 
and  energy  loss  may  only  be  transmitting  from  the  second  row  back  into  the 
first  row. 

In  calculation  906.1020*  the  deposited  energy  was  increased  to  2.364  x  10lc 
ergs  per  cell.  The  transmissive  lower  boundary  was  again  specified.  The 
individual  cell  dimensions  were  retained;  therefore  the  effective  charge 
density  was  2.1611  x  107  ergs/mm3,  which  is  still  lower  than  the  desired  charge 
density.  In  Figure  33  the  maximum  reflected  pressure  of  the  station  next  to 
the  right  boundary  was  25.5  MPa.  This  fact  indicates  an  unexplained  reversal 
in  magnitude  compared  to  the  peak  value  of  49.4  MPa  for  calculation  905.1080 
(Fig.  28)  which  was  run  with  the  smaller  deposited  energy.  Pressure  contours 
and  the  velocity  field  for  the  approximate  time  of  peak  pressure  are  shown  in 
Figure  34.  The  pressure  histories  in  Figure  35  shed  some  light  in  explaining 
this  reversal.  The  grid  used  in  all  three  calculations  was  too  coarse. 
Instabilities  are  present  in  the  pressure  histories  for  calculation  905.1080 
at  Stations  18  and  29  adjacent  to  the  plane  of  symmetry.  They  are  even  more 
evident  at  the  same  stations  for  calculation  906.1020.  It  is  conceivable 
that  the  increased  energy  in  calculation  906.1020  was  enough  to  generate  the 
observed  oscillations  which  may  in  turn  have  caused  excessive  attenuation  of 
the  wave. 

The  computer  costs  for  each  of  these  problems  were  $10  to  define  and 
initialize  the  mesh,  $50  to  perform  the  actual  calculation,  and  $20  to 
plot  the  output.  If  the  mesh  is  refined  the  cost  would  increase  accordingly, 
perhaps  as  much  as  an  order  of  magnitude. 

The  results  of  these  calculations  were  not  fully  conclusive  in  that  they  do 
not  yet  provide  a  basis  for  comparison  with  the  calculation  of  the  1-D 
simultaneous  detonation. 

♦Hydrocode  Calculation  No.  906.1020,  Data  Report,  UNM/CERF  AST-38,  December 
1978. 
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Figure  35.  HULL  Calculations  905.1080  and  906.1020--Instabil ities  in 
Calculations 


SECTION  IV 


DABS  SIMULATOR  CALCULATIONS 

The  DABS  simulator  consists  in  general  of  a  partially  buried  shock  tunnel 
in  which  an  explosive  charge  placed  against  the  closed  end  is  used  to 
produce  a  planar  shock  wave  which  propagates  down  the  tunnel  and  produces  a 
dynamic  flow  behind  the  shock.  Items  of  interest  in  the  calculational  modeling 
of  the  DABS  experiment  include  the  proper  method  for  modeling  the  detonation 
(a  burn  algorithm  versus  an  isothermal  energy  dump),  the  selection  of  appro¬ 
priate  cell  sizes  and  boundary  conditions,  and  the  representation  of  target 
structures  and  reflecting  surfaces. 

A  number  of  1-D  and  2-D  calculations  were  performed  with  the  use  of  the  SAP 
and  HULL  hydrocodes  to  investigate  these  and  other  items  of  interest.  Several 
simulator  configurations  were  then  modeled  and  calculated. 


PARAMETRIC  STUDY  OF  BLAST  WAVE  CHARACTERISTICS 
AS  A  FUNCTION  OF  CELL  SIZE  AND  EXPLOSIVE 
CHARGE  DENSITY  (SAP  HYDROCODE) 


The  results  of  a  series  of  SAP  1-D  numerical  simulations  are  presented  here 
to  provide  a  data  bare  which  can  be  used  to  derive  an  economical  technique 
for  predicting  blast  wave  characteristics  for  DABS  experiments. 

The  fluid  behavior  for  these  simulations  was  governed  by  the  AFWL  variable- 
specific-heat  equation  of  state  for  air.  In  all  cases  Eulerian  flow  was 
simulated  in  a  uniform  mesh  within  a  Cartesian  geometry  and  in  the  absence 
of  artificial  viscous  damping.  The  left  boundary  was  coincident  with  the 
plane  of  symmetry.  Therefore  these  solutions  describe  the  propagation  of 
planar  blast  waves  in  a  closed,  constant-area  tube  which  initially  contains 
air  at  one  standard  atmosphere. 

The  blast  waves  were  generated  by  incrementing  the  internal  energy  in  the  cell 
adjacent  to  the  plane  of  symmetry  between  the  first  and  second  calculational 
cycles.  Consequently  this  first  cell  is  comparable  to  the  driver  section  of  a 
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closed  shock  tube  where  both  it  and  the  reaction  section  contain  air. 

Six  calculations*  were  performed  in  which  both  the  initial  planar  energy 
density  of  the  driver  and  the  cell  size  were  varied.  The  calculational 
matrix  specified  cell  widths  of  10,  100,  and  200  mm  at  energy  levels  of 
10s5  and  1011  ergs/mm2  (Fig.  36).  In  each  calculation  the  histories  of 
hydrodynamic  variables  were  recorded  at  ranges  of  1 ,  2,  3,  and  4  tn  from  the 
plane  of  symmetry.  In  most  instances  it  is  desired  to  predict  the  blast 
intensity  at  some  target  structure.  Therefore  these  data  were  recorded  at 
the  ce1!  immediately  upstream  from  the  nominal  station  range  (Fig.  37)  and 
represent  fluid  properties  at  the  center  of  the  cell. 

When  assessing  structure  response  to  impulsive  blast  loadings,  four  parameters 
are  of  particular  interest: 

1 .  wave  velocity, 

2.  rate  at  which  the  structure  is  loaded, 

3.  maximum  applied  load,  and 

4.  total  applied  load. 

Some  typical  calculated  results  are  presented  in  Figure  38  where  the  waveforms 
for  the  same  location  are  compared  for  three  cell  sizes. 

These  parameters  are  presented  graphically  as  families  of  functions  represent¬ 
ing  appropriate  calculated  blast  wave  characteristics,  viz: 

1.  time  of  arrival  (Fig.  39), 

2.  rise  time  (Fig.  40) , 

3.  peak  pressure  (Figs.  41  and  42),  and 

4.  pressure  impulse  (Fig.  43). 

The  accuracy  with  which  a  nuclear  airblast  environment  can  be  simulated 

is  to  a  great  extent  dependent  on  how  well  these  parameters  can  be  controlled. 


*Hydrocode  Calculations  Nos.  908.0011  through  908.0061,  Data  Report,  UNM/CERF 
AST-44,  December  1978. 
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Figure  38.  Comparison  of  Peak  Pressures  and  Waveforms  for  Different  Cell  Sizes 
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At  ranges  up  to  3  in  shock  arrival  times  exhibited  a  nearly  linear  decay  as 
cell  width  increased  (Fig.  39).  These  solutions  indicate,  however,  that  at 
4  in  this  relation  no  longer  holds,  particularly  for  lower  strength  shocks 
(peak  pressures  <1  MPa)  where  the  effect  of  cell  size  on  shock  velocity  may 
differ  significantly. 

Some  difficulty  was  encountered  in  determining  rise  times  from  pressure 
history  plots.  However,  the  graphic  display  of  these  data  sets  (Fig.  40) 
reveals  a  reasonably  uniform  influence  from  an  increase  of  cell  width.  At  all 
ranges  the  rate  at  which  rise  times  lengthen  gradually  diminishes.  As  one 
might  expect  cell  width  has  the  least  impact  at  closer-in  ranges  where  propa¬ 
gation  velocities  are  greatest. 

The  graphic  interpretation  of  peak  pressure  data  was  based  on  the  presumption 
that  as  the  cell  size  becomes  infinitely  small,  the  calculated  peak  must 
asymptotically  approach  some  limiting  maximum  value.  Therefore,  when  consid¬ 
ering  peak  values  as  a  function  of  cell  width  (Fig.  41),  it  was  necessary  to 
arbitrarily  introduce  an  inflection  in  the  curve  to  satisfy  this  condition.  Due 
to  the  character  of  the  data  in  these  limited  samplings  (three  cell  widths), 
it  is  implied  that  as  the  cell  becomes  very  large  the  peak  pressure  would 
approach  a  lower  limit.  It  is  reasonable  to  contend  that  in  reality  this  is 
the  case.  At  both  shock  strengths  it  is  clear  that  as  range  increases,  peak 
pressure  becomes  markedly  less  dependent  on  cell  width.  The  displays  of  peak 
pressures  versus  range  (Fig.  42)  clearly  show  that  as  cell  width  increases, 
a  rapid  convergence  to  a  minimum  peak  pressure  is  nearly  independent  of 
range. 


When  considering  the  presentations  of  pressure  impulse  (Fig.  43),  one  should 
keep  in  mind  that  the  accumulations  represent  relatively  short  time  intervals, 
35  ms  .and  85  ns  for  energy  levels  of  108  and  10' 1  ergs/mm2,  respectively. 

In  the  extreme  this  means  that  the  pressure  had  decayed  to  only  57  percent  of 
the  peak.  Nevertheless  there  is  an  indication  that  the  impulse  is  influenced 
by  cell  size,  increasing  slightly  at  first  and  then  decreasing  when  the  cell 
width  exceeds  150  mm. 
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Before  using  these  results  to  assess  the  effects  of  cell  size  on  blast  wave 
characteristics,  one  should  consider  certain  aspects  of  this  calculational 
study. 


First,  the  simplistic  technique  used  to  generate  the  airblasts  should  be 
evaluated  with  respect  to  the  simulation  of  the  detonation  of  an  explosive. 

It  is  quite  possible  that  a  multimaterial  model  which  incorporates  the 
effects  of  the  detonation  process  could  significantly  alter  these  solutions. 

Second,  the  graphical  interpretation  of  functions  based  on  data  representing 
only  three  cell  widths  must  be  considered  questionable  in  some  respects.  For 
example,  to  verify  the  assumptions  made  in  regard  to  peak  pressure  behavior, 
a  broader  statistical  base  is  required.  Finally,  since  DABS  tests  may  involve 
tunnel  lengths  well  in  excess  of  10  m,  calculations  of  sufficient  duration 
to  obtain  data  at  appropriate  ranges  are  recommended. 

SIMULATION  OF  4.14-MPa  REFLECTED  AIRBLAST  WAVE 

SAP  calculation  906.1081*  was  performed  to  evaluate  the  incident  airblast  wave 
in  an  explosive  simulation  of  a  1 25- kt  nuclear  airblast  wave  in  support  of 
the  early  HAVE  HOST  T  series  (in-trench)  tests.  A  range  was  sought  which 
produces  a  peak  pressure  of  4.14  Mpa.  Calculation  906.1091**  then  evaluated 
the  flow  with  a  solid,  reflecting  boundary  introduced  into  the  shock  tunnel 
at  the  4.14-Mpa  range  found  in  calculation  906.1081. 

The  driver  consisted  of  25  cells  of  burned  ammonium  nitrate  and  fuel  oil  (ANF0) 
with  a  mass  density  of  1.38  x  10'3  g/mm3  and  an  energy  density  of  2.822  x  1010 
ergs/g  as  shown  in  Figure  44.  Cell  thickness  was  approximately  3.6232  mm  for 
a  total  driver  width  of  90.5797  mm.  The  total  driver  mass  was  12.5  g,  and 
total  driver  energy  was  3.5275  x  1011  ergs.  The  shock  tunnel  contained  real 
air  with  ambient  densities  of  1.225  x  10”  g/mm3  and  2.04448  x  10'  ergs/g 
and  a  cell  thickness  of  250  mm.  SAP  method  Number  1  (Lagrangian)  was 
selected  as  the  hydrocode  to  be  used  because  of  its  multimaterial  capability 


*Hydrocode  Calculation  No.  906.1081,  Data  Report,  UNM/CERF  AST-26,  December 
1  978. 

**Hydrocode  Calculation  No.  906.1091,  Data  Report,  UNM/CERF  AST-27,  December 
1  978. 
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Incident  Wave  Geometry--Calculation  906.1081 
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Reflected  Wave  Geometry--Calculation  906.1091 
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Figure  44.  Geometrical  Descriptions  for  Calculations  906.1081  and  906.1091 


and  because  its  Lagrangian  properties  would  allow  the  interaction  of  detonation 
products  to  be  traced  throughout  the  calculations.  The  displacement  of  this 
interface  at  various  times  for  both  incident  and  reflected  wave  calculations 
is  recorded  in  Table  4.  An  artificial  viscosity  function  (Ref.  12)  was 
employed  to  provide  a  stable,  sharply  defined  shock  front.  The  mesh  geometries 
for  the  two  calculations  are  identical  in  the  driver  region  and  also  in  the 
reaction  region  of  the  shock  tunnel.  Differences  occur  at  the  reflecting 
boundary  which  is  located  at  65,500  mm  in  calculation  906.1091  as  shown  in 
Figure  44.  Material  and  energy  deposition  are  identical  for  the  two  calcula¬ 
tions  . 

Pressure  history  stations  were  set  in  calculation  906.1081  at  1-m  intervals, 
starting  at  57  m  and  continuing  through  79  m.  The  calculation  was  then  run  far 
enough  to  identify  that  the  peak  pressure  reached  a  value  of  4.14  MPa  at  the 
65-m  range  at  approximately  25.5  ms  (Figs.  45  and  46).  It  was  then  necessary 
to  extend  the  open  end  of  the  tunnel  and  allow  the  calculation  to  continue 
until  the  interface  of  detonation  products  passed  well  beyond  the  station  of 
interest.  Passage  of  this  interface  may  be  noted  in  Figure  45  as  the  change 
in  slope  of  the  pressure  curve  at  approximately  45.5  m.  It  is  also  verified 
in  Figure  47  as  the  abrupt  change  in  density  at  65  m.  The  pressure  history 
stations  for  calculation  906.1091  were  complemented  with  additional  stations 
at  100-mm  intervals  starting  at  64.5  m  through  65.4  m  to  measure  the  reflected 
pressure.  Calculation  906.1091  was  run  to  60  ms  to  allow  sufficient  time  fur 
the  detonation  products  to  interact  with  the  reflected  wave.  The  reflected 
wave,  traveling  to  the  left,  intercepts  the  denser  detonation  products  which 
are  still  moving  to  the  right,  resulting  in  a  secondary  reflection  wave  toward 
the  right.  The  moment  of  interception  can  be  seen  in  Figure  48  in  the  spikes 
at  a  radius  of  56  m.  This  secondary  reflection  is  then  reflected  in  turn  at 
the  cavity  boundary  and  appears  as  the  second  peak  of  pressure  in  Figure  49 
with  a  maximum  value  of  16  MPa  at  approximately  44  ms.  The  moment  of  arrival 
of  the  secondary  wave  at  65.4  m  is  shown  in  Figure  50.  The  interface  boundary 
reaches  a  maximum  radius  of  56.24  m  at  39  ms  (Table  4)  before  it  starts  con¬ 
tracting  after  the  passage  of  the  original  reflected  wave. 
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TABLE  4.  POSITION  OF  DETONATION  PRODUCTS  BOUNDARY 


Time,  ms 

Incident  Uave 
Radius,  m 

Reflected  Wave 
Radius,  m 

0 

0.0906 

0.0906 

5 

13.6949 

13.6949 

10 

23.5236 

23.5236 

15 

31 .6707 

31.6707 

20 

38.7236 

38.7236 

25 

44.9707 

44.9707 

30 

50.5805 

50.5790 

35 

55.6693 

55.6623 

40 

60.3230 

56.2114 

45 

64.5855 

55.5395 

50 

68.5237 

54.0469 

55 

72.1597 

51.1936 

60 

74.2061 

48.8432 

Each  calculation  was  run  to  a  problem  time  of  60  ms.  Calculation  906.1081  re 
quired  5125  cycles  and  cost  $75.  Calculation  906.1091  required  9723  cycles 
and  cost  $120. 

The  presence  of  the  second  reflected  peak  has  been  noted  in  previous  APOD 
calculations  used  in  predicting  the  test  environments  for  the  HAVE  HOST  T-l 
Event.*  The  corresponding  experimental  measurements  of  reflected  pressure 
at  the  plug  of  the  T-l  event**  only  hints  at  the  presence  of  this  second 
reflection.  It  should  be  noted  that  the  interface  in  the  calculations  is  a 
sharply  defined  boundary  whereas  in  reality  the  interface  becomes  diffused 
and  probably  does  not  stand  out  to  the  extent  noted  in  the  calculations. 

This  conclusion  was  confirmed  by  an  Eulerian  HULL  calculation  (reported  later 
in  the  section). 


*Renick,  J. ,  et  al . ,  HAVE  HOJT  T-l  Predictions,  Letter  Report,  Civil 
Engineering  Research  Division,  Air  Force  Weapons  Laboratory,  Kirtland  Air 
Force  Base,  New  Mexico,  12  July  1977. 

**UAVii  liodT  T-l  >,:<ick  Look,  Preliminary  Report,  Civil  Engineering  Research 
Division,  Air  Force  Weapons  Laboratory,  Kirtland  Air  Force  Base,  New  Mexico, 
25  August  1 977 . 
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Calculation:  906.1081 
Radius:  65  m 

Maximum:  4.1683  MPa 
Time  Of  Maximum:  25.513 


Time,  ms 
ressure  Impulse 


Calculation:  906.1081 
Radius:  65  m 
Maximum:  0.049827  MPa-s 
Time  of  Maximum:  59.675 


Density 


Fiqure  46.  Incident  Wave  Profiles  at  26  ms 
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Figure  47.  Incident  Wave  Profiles  at  46  ms 
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PARAMETRIC  STUDY  OF  SAP  LAGRANGIAN  HYDROCODE 

In  the  course  of  varying  the  mass  and  energy  deposition  to  achieve  the 
4.14-MPa  incident  pressure  loading  during  the  previous  calculation  for 
the  HAVE  HOST  T  series,  some  unexpected  results  having  to  do  with  the  number 
of  cells  used  to  define  the  explosive  driver  section  and  the  mass  and  energy 
in  each  were  noted.  Not  all  of  the  results  are  desirable,  and  it  is  important 
to  be  able  to  predict  their  occurrence.  Seven  abbreviated  trial  calculations* 
were  made  in  the  course  of  finding  the  proper  mass-energy  deposition  for  the 
driver  section  to  achieve  the  4.14-MPa  wave.  Variations  were  made  in  1)  the 
number  of  cell  divisions  within  the  driver,  2)  the  total  width  and  mass  density 
of  the  driver,  and  3)  the  width  of  the  air  cells  in  the  shock  tunnel.  Table  5 
provides  a  listing  of  all  the  input  parameters  as  well  as  a  comparative  look 
at  the  peak  pressure  at  the  5-m  station.  Three  distinct  energy  levels  were 
used  in  the  study  with  three  calculations  at  each  level.  Common  to  all  calcula¬ 
tions  were  the  use  of  the  real  air  EOS  in  the  reaction  cavity,  the  same  multi¬ 
material  Lagrangian  version  of  the  SAP  hydrocode,  and  the  burned  ANFO  EOS  in 
the  driver  section.  The  ANFO  EOS  did  have  modified  mass  and  energy  densities 
as  listed  in  Table  5. 

Because  previous  calculations  had  employed  an  energy  dump  into  a  single 
driver  cell,  this  procedure  was  followed  in  the  initial  attempts  in  this  set 
of  calculations.  It  was  noted,  however,  that  the  dump  into  the  single  cell 
resulted  in  unstable  oscillations  during  the  late  times  after  the  arrival  of 
the  initial  peak  and  during  the  subsequent  pressure  decay.  These  oscillations 
were  noted  in  calculations  906.1011  (Fig.  51),  906.1021,  and  906.1031.  Another 
observation  made  from  the  use  of  the  single-driver  cell  was  the  lack  of  full 
downstream  dissipation  of  the  driver  energy.  When  multiple  cells  were  em¬ 
ployed  the  peak  pressure  at  downstream  stations  seemed  to  peak  when 


*Hydrocode  Calculations  Nos.  906.1011  through  906.1071  Data  Report,  UN'M/CERF 
AST-30,  December  1978. 
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TABLE  5.  EFFECTS  OF  CALCULATION  PARAMETERS  ON  PEAK  PRESSURE 
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approximately  10  cells  were  used  in  the  driver.  (Compare  calculations  906.1011 
versus  906.1041  in  Table  5  and  then  906.1021  versus  906.1061.) 

It  was  also  observed  that  cell  size  within  the  reaction  cavity  had  a  definite 
effect  on  the  magnitude  of  the  peak  pressure  produced  downstream.  Figures 
52,  53,  and  54  are  plots  of  peak  pressure  versus  station  radii  for  the  three 
groups  of  energy  levels.  The  energy  level  is  constant  for  all  calculations 
within  a  group.  The  legend  on  the  figures  identifies  the  particular  curve 
according  to  its  pressure  at  the  5-m  location.  In  all  three  group  plots 
the  highest  peak  pressure  occurs  in  the  calculation  with  the  smallest  air 
cell  width  and  decreases  to  the  lowest  peak  with  the  largest  cell  width. 

Of  special  interest  are  calculations  906.1081  and  906.5081.*  Calculation 
906.1081  is  the  final  version  used  in  the  4.14-MPa  wave  calculation.  Calcula¬ 
tion  906.5081  uses  the  identical  input  parameters  as  calculation  906.1081  with 
the  sole  exception  of  using  a  burn  subroutine  to  release  the  driver  energy 
rather  than  using  an  energy  dump.  In  Figure  53  the  peak  pressures  at  all 
stations  for  calculation  906.5081  are  from  30  to  65  percent  higher  than  the 
corresponding  stations  for  calculation  906.1081.  In  both  the  energy  dump  and 
burn  methods  a  pressure  differential  exists  initially  between  the  last  driver 
cell  and  the  first  air  cell.  The  burn  subroutine  begins  adding  energy  at  the 
other  end  of  the  driver  which  establishes  a  second  pressure  front  at  the  rear 
of  the  driver  section.  This  leads  to  the  possibility  of  an  early  reflection 
off  the  backwall  to  boost  the  outgoing  pressure  wave  above  the  values  noted 
in  the  no-burn  process. 

From  these  parametric  calculations  it  can  be  concluded  that  a  simple  iso¬ 
thermal  energy  dump  is  not  adequate  for  the  accurate  modeling  of  a  detonation 
process  occurring  across  some  significant  thickness  of  explosive  material.  It 
can,  however,  be  used  for  qualitative  calculations  where  an  exact  modeling  of 
the  energy  source  is  not  important. 

*Hydrocode  Calculation  No.  906.1081,  Data  Report,  UNM/CERF  AST-26, 

December  1978. 

Hydrocode  Calculation  No.  906.5081,  Data  Report,  UNM/CERF  AST-28, 

December  1978. 
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Figure  52.  Peak  Pressure  versus  Length  (G 
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Figure  53.  Peak  Pressure  versus  Length  (Group  II) 


Calculation:  906.1021 
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Figure  54.  Peak  Pressure  versus  Length  (Group  III) 


EVALUATION  OF  DIFFUSION  LIMITATION  IN  HULL  CODE 


Airblast  waves  generated  in  explosively  driven  shock  tunnels  and  test  cavities 
are  essentially  high-strength  planar  shock  waves  which,  it  can  be  assumed, 
are  representative  of  the  interior  ballistics  encountered  in  typical  shock 
tube  experiments.  The  contact  surfaces*  which  accompany  shocks  of  this  kind 
can  play  an  important  part  in  complex  shock  interactions  during  the  evolution 
of  the  wave  system.  The  purpose  of  the  numerical  simulation  described  here 
(calculation  906.0040**)  was  to  evaluate  the  diffusion  limitation  scheme  which 
is  automatically  implemented  in  HULL  for  multimaterial  models.  Since  the 
HULL  code  simulates  flow  in  a  Eulerian  frame  of  reference,  some  mechanism 
of  this  kind  is  required  if  contact  discontinuities  are  to  be  realistically 
preserved. 

A  model  was  chosen  which  was  readily  amenable  to  a  closed  form  analytical 
solution.  In  this  instance  a  20-m  closed  shock  tube  with  the  diaphragm 
at  midlength  was  simulated  (Fig.  55).  With  these  proportions  the  primary 
shock  and  contact  surface  can  diverge  suitably  before  any  reflections  take 
place  (Fig.  56).  This  configuration  was  modeled  in  a  Cartesian  coordinate 
system  in  which  the  closed  end  of  the  high  pressure  chamber  was  at  the 
plane  of  symmetry.  Since  flow  in  this  case  is  one-dimensional,  only  the 
minimum  number  of  cells  (6)  required  by  the  HULL  code  were  used  in  the  trans¬ 
verse  or  "radial"  direction.  All  cells  were  100  by  100  mm,  and  all  boundaries 
were  reflective. 

Helium,  pressurized  to  100  atmospheres,  was  selected  as  the  driver  gas  and 
air  at  standard  sea  level  conditions  as  the  reaction  gas.  Initial  conditions 
for  the  driver  (Region  4,  Fig.  56)  were  assumed  to  be  the  result  of  an  isen- 
tropic  compression  (Table  6,  Ref.  16).  Instantaneous  removal  of  the  diaphragm 


16.  Liepmann,  H.  W.,  and  Roshko,  A.,  Elements  of  Gas  Dynamics,  pp.  62-83, 
Wiley,  1957. 

*  A  surface  which  separates  the  driver  gas  and  reaction  gas  in  a  shock  tube. 
It  is  in  reality  a  region  but  is  often  theoretically  idealized  to  be  a 
discrete  material  interface. 

**Hydrocode  Calculation  No.  906.0040,  Data  Report,  UNM/CERF  AST-37,  December 
1978. 


102 


Regions  (3)  and  (4)  Contain  Helium 
Regions  (1)  and  (2)  Contain  Air 
Diaphram  Pressure  Ratio  (Pu/Pi)  =  100 
Shock  Strength  (P2/P1)  =  33.66 


Figure  56.  Wave  System  Diaqrain--Calculation  906.0040 
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TABLE  6.  DRIVER  REGION  INITIAL  CONDITIONS3  — 
HULL  CALCULATION  906.0040 


P-  =  Po  (^-)Y  1  =  2.679  x  10 
P„  =  100  P0  =10  MPa 


^4 


[p4  (Y  -  1)] 


5.598  x 


where  P0  =  0.1  MPa 

5 

f  ~  3 

«.:•(  =  1 .6905  x  1  O’7  g/mm? 
Tc  =  2.93c  K 


6  g/mm3 


1010  ergs/g 


K 


Helium  isentropically  compressed  to  100  atmospheres 


separating  Regions  4  and  1  was  simulated  by  using  these  values  to  describe  the 
initial  state  of  the  driver  gas.  A  solution  for  the  fluid  response  of  this 
model  was  recorded  every  100  ps  for  a  period  of  3  ms. 

To  obtain  the  closed  form  analytical  solution  based  on  shock  tube  theory 
the  computer  program  SH0K2B  (Appendix  C)  was  written.  Given  the  initial 
conditions  for  the  driver  and  reaction  gases  (Regions  4  and  1)  and  the  desired 
solution  time,  this  program  first  solves  the  basic  shock  tube  equation  for 
the  shock  strength  (P2/Pi)  which  is  defined  implicitly  as  a  function  of  the 
diaphragm  pressure  ratio  (P4/P1).  The  remaining  state  variables  are  then  cal¬ 
culated  for  Regions  2  and  3.  Finally  the  positions  of  primary  shock  front, 
contact  surface,  and  expansion  wave  are  determined. 

The  HULL  numerical  solution  density,  specific  internal  energy,  and  pressure 
profiles  at  3  ms  is  shown  in  Figures  57  through  59.  For  the  purpose  of 
comparison  the  analytical  solution  has  been  overlayed.  It  can  be  seen  that 
these  solutions  are  in  remarkably  close  agreement.  The  gradients  in  Region  2 
in  both  the  density  and  energy  profiles  (Figs.  57  and  58)  show  the  extent 
to  which  HULL  has  permitted  diffusion  to  take  place  at  the  contact  surface. 

The  HULL  solution  locates  the  contact  surface  approximately  14.1  m  as  opposed 
to  the  analytical  value  of  14.4  m. 

From  this  calculation  it  can  be  concluded  that  the  HULL  diffusion  limitation 
scheme  can  creditably  simulate  the  contact  discontinuity  of  the  simplistic 
shock  tube  environment.  Since  in  reality  such  a  discontinuity  cannot  be 
sustained  (Ref.  17),  a  diffusion  does  occur.  The  truth  then  may  very  well 
lie  somewhere  between  this  HULL  prediction  and  the  theoretical  case  which 
presumes  no  diffusion. 

Although  blast  effects  in  explosively  driven  test  facilities  are  essentially 
analogous  to  the  1-D  wave  motions  and  phenomena  of  the  shock  tube,  certain 
physical  factors  were  overlooked  here.  For  example,  the  explosive  energies 
involved  are  equivalent  to  diaphragm  pressure  ratios  that  may  be  orders  of 
magnitude  greater  than  can  be  achieved  in  a  diaphragm  shock  tube  experiment. 

17.  Bradley,  J.  N.,  Shook  Waves  in  Chemistry  and  Physics ,  pp.  103-4,  Wiley, 
1962. 
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Figure  58.  Axial  Energy  Profile--Calculation  906.0040  at  3  ms 
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Figure  59.  Axial  Pressure  Profile--Calculation  906.0040  at  3  ms 


Also,  the  typical  test  configuration  is  such  that  the  driver  section  is 
relatively  short.  This  condition  will  most  certainly  lead  to  complex  wave 
interactions  involving  reflections,  refractions,  overtaking  situations,  and 
possibly  the  generation  of  contact  regions  of  variable  entropy.  Finally,  of 
course,  it  is  to  be  expected  that  the  explosive  burning  processes  would 
contribute  significantly  to  the  experimental  environment. 

A  number  of  closed-form  analytical  solutions  for  various  wave  interactions 
are  available  (Ref.  18).  It  is  therefore  recommended  that  HULL  simulations  of 
some  of  the  more  straightforward  cases  be  evaluated. 


HULL  NUMERICAL  SIMULATION  OF  NORMAL  REFLECTION 
OF  4. 14- MPa  PLANAR  SHOCK  WAVE 

HULL  calculation  906.1090*  predicts  the  blast  effects  in  an  explosively 
driven  shock  tunnel.  It  was  part  of  a  calculational  study  to  determine  driver 
specifications  for  HAVE  HOST  experiments.  The  specifications  for  this  model 
duplicated  those  for  a  1-D  SAP  Lagrangian  calculation  (calculation  906.1091**). 
For  all  practical  purposes  it  simulates  a  65.4-m  closed  shock  tunnel  which 
is  represented  in  a  Cartesian  coordinate  system  (Fig.  60).  All  boundaries 
of  the  rectangular  mesh  are  reflective. 

The  explosive  driver  initially  occupied  the  column  of  cells  adjacent  to  the 
left  boundary  so  as  to  generate  the  planar  blast  wave.  This  explosive  was 
then  considered  to  be  motionless  and  in  a  completely  burned  state  with  a 
specific  internal  energy  of  2.822  x  1010  ergs/g  and  a  mass  density  of 
8.5  x  10”4  g/mm2.  To  obtain  an  incident  peak  pressure  of  4.14  MPa  at  65.4  m 

4 

the  thickness  of  the  driver  was  set  at  90.5797  mm  (i.e.,  a  planar  energy 
density  of  2.17  x  109  ergs/mm2),  The  behavior  of  this  burned  explosive 
gas  governed  by  the  HULL  EOS  for  burned  ANFO.  The  remainder  of  the  mesh  was 


18.  Glass,  I.  I.,  and  Hall,  J.  G.,  Handbook  of  Supersonic  Aerodynamics ,  Vol .  6, 
Section  18,  NAVORD  Report  1488,  pp.  81-101,  December  1959. 

*  Hydrocode  Calculation  No.  906.1090,  Data  Report,  UNM/CERF  AST-36,  December 
1978. 

**  Hydrocode  Calculation  No.  906.1091,  Data  Report,  UNM/CERF  AST-27,  December 
1978. 
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Boundary  Conditions--Ref1ective 


Initial  Conditions 


Equation  of  State 

AX,  mm 

Internal  Energy, 
ergs/g 

Density,  g/mm3 

Burned  ANFO 

0  to  90.5797 

mansBi 

8.5  x  10'4 

Real  Air 

90.5797  to  65,400 

2.0448  x  109 

1 .225  x  10‘6 

Figure  60.  Model  Specifications--HOLL  Simulation  906.1090 
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defined  as  real  air  (i.e.,  variable  specific  heats)  at  conditions  correspond¬ 
ing  to  1  standard  atmosphere  at  sea  level.  Calculated  peak  pressure  as  a 
function  of  range  is  shown  in  Figure  61. 

This  model  represents  a  typical  case  of  interior  shock  physics  where  the 
exploded  gas  is  separated  from  the  air  in  the  tunnel  by  a  membrane  which  is 
instantaneously  removed  (Ref.  14,  pp.  173-176).  Since  the  contact  surface 
(membrane)  is  quite  close  to  the  closed  end,  the  rarefaction  reflects  quite 
early  and  a  complex  sequence  of  wave  interactions  and  reflections  takes  place. 
Of  particular  importance  here  is  the  overtaking  of  the  primary  shock  by  the 
reflected  rarefaction.  The  outcome  of  such  an  interaction  is  dependent  on  the 
relative  strength  of  the  rarefaction.  Glass  and  Hall  (Ref.  18)  have  addressed 
this  subject  and  described  the  various  wave  systems  which  can  result.  In  all 
cases,  the  attenuation  of  the  shock  gives  rise  to  a  variable  entropy  field 
in  the  form  of  a  contact  region  rather  than  a  contact  surface. 

Pressure  histories  at  60  and  65.4  m  (Figs.  62  and  63)  depict  the  character 
of  the  incident  and  reflected  shock  waves.  A  peak  reflected  pressure  at 
the  right  boundary  of  39  MPa  is  equivalent  to  a  peak  pressure  ratio  of  9.4. 
This  value  is  considerably  higher  than  the  theoretical  (7.5);  however, 
it  should  be  noted  that  this  calculation  did  not  include  any  artificial  vis¬ 
cous  damping.  The  pressure  profile  just  prior  to  incidence  (Fig.  64)  clearly 
shows  the  erosive  effects  of  overtaking  rarefactions.  This  waveform  can  be 
compared  directly  with  Figure  46  which  clearly  reveals  the  presence  of  the 
original  contact  surface  which  is  rigorously  preserved  by  the  Lagrangian 
calculation.  Experimental  data  show  this  Eulerian  HULL  calculation  to  be 
more  representative  of  real  world  effects. 

The  HULL  code  diffusion  limiter  would  tend  to  preserve  the  original  contact 
surface.  However,  since  the  driver  section  is  very  short  there  is  ample  time 
for  the  development  of  a  complex  wave  system  involving  a  succession  of  re¬ 
flections  and  refractions  as  one  wave  overtakes  another.  These  interactions 
can  lead  to  the  creation  of  additional  discontinuities.  It  is  somewhat 
less  than  clear  as  to  the  extent  this  solution  reflects  either  theory  or 
rea 1 i ty . 
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Figure  62.  Pressure  History  at  60  m--Calculation  906.1090 
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Any  evaluation  must  take  into  account  another  aspect  of  blast  environments 
encountered  in  explosively  driven  shock  tunnels.  In  this  simulation,  for 
example,  the  initial  pressure  in  the  driver  section  is  approximately  6000  MPa. 
This  is  equivalent  to  a  diaphragm  pressure  ratio  of  60,000.  It  is  not  at  all 
certain  that  classical  shock  tube  theory  is  valid  in  these  higher  pressure 
regimes  and  it  may  well  be  that  numerical  solutions  such  as  this  can  only  be 
analyzed  in  the  light  of  the  results  from  well -instrumented  and  controlled 
high-pressure  experiments. 

If  one  makes  a  reasonable  allowance  for  the  absence  of  artificial  viscous 
damping,  an  estimated  reflected  peak  pressure  of  30  MPa  compares  favorably 
with  a  analytical  value  of  31.05.  Finally,  examination  of  the  pressure, 
density,  and  energy  profiles  representing  the  solution  at  20  ms  (Figs.  64, 

65,  and  66)  reveals  a  discontinuity  at  40  m.  This  may  indicate  the  presence 
of  a  contact  region  which  is  generated  as  a  rarefaction  erodes  the  shock 
wave. 


HULL  SIMULATION  OF  TYPICAL  DABS  EXPERIMENT 

HULL  calculation  908.0010*  was  an  attempt  to  simulate  the  blast  effects 
numerically  which  are  typical  of  those  experienced  during  recently  fielded 
large  DABS  experiments.  The  DABS  test  facility  is  an  explosively  driven, 
partially  buried  shock  tunnel  which  is  used  to  impose  dynamic  airblast 
loadings  on  scaled  target  structures.  This  calculation  was  part  of  a  parame¬ 
ter  study  to  establish  design  guidelines  to  be  used  in  preparing  specifica¬ 
tions  for  future  experiments. 

The  working  medium  for  this  simulation  was  modeled  as  air,  the  behavior  of 
which  was  governed  by  the  AFWL  variable  specific  heat  EOS  (Ref.  19).  The 
shock  tunnel  was  represented  in  a  Cartesian  coordinate  system  as  the  space 
between  two  parallel  and  reflective  X-Z  planes  which  are  4.1  m  apart  (Fig.  67). 
Consequently  this  solution  simulates  only  a  2-D  fluid  response. 


19.  Needham,  C.  E.,  Nualcai •  Blast  Btandai'd  (l  kt) ,  AFWL-TR-73-55  (Rev.),  Air 
Force  Weapons  Laboratory,  Kirtland  Air  Force  Base,  New  Mexico,  1975. 

*  Hydrocode  Calculation  No.  908.0010,  Data  Report,  UNM/CERF  AST-42,  December 
1978. 
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Figure  65.  Axial  Density  Profile  at  20  ms--Calculation  906.1090 


118 


Energy  erqs/q  (x  10  M 


The  explosive  driver  was  located  at  the  closed  end  of  the  tunnel  (plane  of 
symmetry)  and  was  initially  configured  as  being  200  mm  thick  from  floor  to 
ceiling.  The  open  end  of  the  tunnel  is  simulated  with  a  transmissive  boundary 
which  is  16  m  from  the  plane  of  symmetry.  The  rectangular  target  structure 
is  modeled  by  a  1-  by  6-m  "island"  (a  region  of  the  mesh  which  is  filled  with 
a  motionless  fluid  of  infinite  strength).  In  short  all  tunnel  and  structure 
surfaces  are  modeled  as  perfectly  rigid  material  interfaces  (i.e.,  perfectly 
reflecti ve) . 

The  instantaneous  release  of  explosive  energy  was  simulated  by  a  one-time 
deposition  between  the  second  and  third  calculational  cycles.  This  energy 
deposition  was  made  in  the  column  of  cells  adjacent  to  the  plane  of  symmetry 
and  corresponds  to  a  planar  energy  density  of  1.195  x  1 0 n  ergs/mm2. 

The  energy  deposition  was  chosen  to  generate  a  blast  wave  which  would  deliver 
a  peak  incident  pressure  of  approximately  3.448  MPa  at  a  range  corresponding 
to  the  vertical  front  face  of  the  target  structure  (10  m).  Dimensions  for 
the  tunnel  as  well  as  the  position  of  the  structure  were  patterned  after  the 
HAVE  HOST  S-l  experiment. 

This  model  is  essentially  analogous  to  the  case  for  a  nonuniform  (decaying) 
shock  where  the  front  is  overtaken  by  a  rarefaction  which  remains  in  contact 
indefinitely.  When  the  instantaneous  energy  deposition  is  made,  the  closed 
end  of  the  tunnel  (plane  of  symmetry)  for  all  practical  purposes  represents  a 
piston  which  has  suddenly  stopped.  This,  of  course,  gives  rise  to  an  expansion 
wave  which  at  0.2  ms  (the  approximate  shock  arrival  time  at  2  m)  has  already 
overtaken  and  begun  to  erode  the  primary  shock.  By  the  time  the  shock  has 
reached  the  target  structure  (10  m)  the  peak  incident  pressure  has  decayed 
to  3.73  MPa.  This  calculated  rate  of  decay  compares  quite  well  with  the 
analytical  solution  (Fig.  68)  where  the  shock  strength  is  inversely  propor¬ 
tional  to  the  square  root  of  time. 

The  shock  wave  impinges  on  the  structure  front  face  at  normal  incidence  and 
results  in  a  calculated  peak  reflected  pressure  at  the  base  (Station  25)  of 
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Figure  68.  Incident  Pressure  versus  Range 


23.1  MPa.  This  corresponds  to  a  peak  pressure  ratio  of  6.2  as  opposed  to 
a  more  analytically  exact  value  of  7.5  (Fig.  69).*  This  discrepancy  is 
believed  due  principally  to  the  unrealistic  prescription  of  properties  behind 
the  shock  front.  The  arbitrary  deposition  of  energy  can  give  rise  to  momen¬ 
tum  and  thermodynamic  states  which  are  not  appropriate  for  shocks  in  real  air. 

Diffraction  of  the  planar  blast  wave  is  accompanied  by  a  rarefaction  which 
propagates  downward  from  the  top  corner  of  the  front  face  of  the  target 
structure.  As  a  result  reflected  peak  pressures  on  this  surface  are  signifi¬ 
cantly  attenuated.  For  example,  the  calculated  peak  at  the  corner  (Station 
27)  was  13  MPa  or  nearly  43  percent  less  than  the  value  recorded  at  floor 
level.  At  a  height  of  450  mm  (Station  26)  a  peak  of  21.6  MPa  was  observed. 
This  gradient  is  due  in  large  measure  to  the  pseudoviscous  damping  of  the 
relatively  course  mesh  (200-  by  100-mm  cells).  Rise  times  of  approximately 
0.5  ms  permit  significant  erosion  of  reflections  at  this  surface  (Fig.  70). 

As  the  blast  continues  to  engulf  the  structure,  heated  fluid  from  the  high- 
pressure  region  in  front  of  the  structure  flows  upward  to  develop  a  secondary 
oblique  shock  wave  which  impinges  on  the  tunnel  roof  at  approximately  6.2  ms 
(Fig.  71a).  The  resulting  reflected  shock  arrives  at  the  top  surface  of  the 
structure  at  approximately  8.5  ms  (Fig.  71b).  Peak  reflected  pressures  at 
this  surface  range  fro"'  71  to  82  percent  of  the  incident  peaks  (Fig. 72). 

The  simplistic  model  used  for  this  calculation  precludes  all  but  a  qualitative 
evaluation  of  results.  Some  real-world  conditions  and  phenomena  that  could 
significantly  affect  the  solution  which  were  not  considered  are 

1.  Explosive  materials  and  processes 

2.  Soils  and  structure  strains 

3.  Three-dimensional  (side)  effects 

Of  particular  interest  in  DABS  tests  is  the  likelihood  and  extent  of  target 
structure  loadings  which  are  due  solely  to  the  proximity  of  interior  tunnel 

*  Chown,  W.  H.,  "Calculation  of  Regular  Reflection  of  Oblique  Shocks  Using 
the  AFWL  Equation  of  State  for  Air,"  Internal  Memorandum  to  Lt.  Col.  J.  J. 
Osborn,  Air  Force  Weapons  Laboratory,  Kirtland  Air  Force  Base,  New  Mexico, 
September  7,  1973. 
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Figure  71.  HULL  Simulation— DABS  Parameter  Study 
(Calculation  908.0010)  (2  of  2) 
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Pressures  at  Top  of  Target  Structure 


surfaces.  On  the  presumption  that  in  actuality  dimensional  changes  during  the 
diffraction  phase  are  minimal  this  solution  indicates  that  reflected  shocks 
within  such  a  test  facility  may  contribute  substantially  to  the  airblast 
environment.  Reflections  from  the  tunnel  ceiling  have  caused  repressurizations 
at  the  top  of  the  target  structure  which  correspond  to  impulse  load  increases 
in  the  neighborhood  of  20  percent. 

The  low  peak  reflected  pressures  obtained  with  this  numerical  simulation 
suggest  that  other  characteristics  of  the  solution  may  not  provide  a  reliable 
data  base  for  establishing  design  parameters.  It  is  recommended  that  the 
model  specifications  be  revised  to  more  realistically  represent  the  experi¬ 
mental  environment  insofar  as  is  practical.  In  addition  to  refining  the 
material  models  and  soil-air  interface  behavior  a  significantly  finer  mesh 
may  be  necessary. 
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SECTION  V 

CONCLUSIONS  AND  RECOMMENDATIONS 


The  state  of  the  art  in  modeling  explosive  simulation  experiments  with 
hydrodynamic  computer  codes  has  advanced  markedly  during  the  past  several 
years.  Several  calculational  tools  are  available  ranging  from  simple  1-D 
similarity  solutions  to  complex  2-D  Eulerian  hydrocodes  with  a  variety  of 
geometries  and  boundary  conditions. 

Methods  for  modeling  the  HEST  simulator  have  increased  the  qualitative  under¬ 
standing  of  the  complex  explosion  mechanisms  and  shock  interactions  which 
occur  within  the  HEST  explosion  cavity.  HEST  designers  can  now  use  these 
calculational  tools  to  better  understand  where  to  place  explosive  arrays  and 
instrumentation  for  improved  HEST  performance.  For  example,  the  most  uniform 
and  widely  dispersed  arrangement  of  the  explosive  material  that  is  practical 
should  be  used  when  designing  a  HEST  explosion  cavity.  Instruments  should  be 
placed  at  locations  which  avoid  the  shock  focusing  effects  from  individual 
strands  of  detonating  cord. 

Quantitative  HEST  prediction  models  will  require  the  use  of  an  equation  of 
state  for  the  plastic  foam  materials  used  to  form  the  HEST  explosion  cavity 
(which  is  currently  under  development)  and  a  more  precise  means  for  modeling 
the  release  of  explosive  energy  into  the  HEST  explosion  cavity.  "Burn"  models 
are  not  currently  available  for  either  detonating  cord  or  for  ammonium 
nitrate-based  "slurry"  explosives.  The  use  of  an  isothermal  energy  deposition 
has  proven  to  be  quantitatively  inadequate  as  a  means  of  modeling  explosion 
processes.  It  can,  however,  be  used  for  the  qualitative  evaluation  of  such 
parameters  as  geometrical  effects,  shock  reflections  from  various  structures, 
boundary  conditions,  and  so  on. 

The  DABS  simulator  can  be  realistically  modeled  using  the  above  mentioned 
calculational  tools.  Again,  as  in  HEST,  quantitative  predictions  must  await 
the  development  of  accurate  "burn"  routines  for  the  nonideal ized  explosive 
materials  used  for  the  DABS  driver.  Qualitative  calculations,  however,  are 
very  useful  for  describing  the  shock  interaction  with  downstream  target 
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structures.  Calculated  peak  pressures  and  shock  rise  times  are  quite  dependent 
on  the  cell  sizes  used.  Increasing  the  number  of  cells,  however,  substantially 
increases  the  cost  of  the  calculation.  The  prudent  selection  of  cell  size 
and/or  mesh  rezoning  for  a  proposed  calculation  is  essential. 

Computer  costs  for  hydrocode  calculations  of  explosive  simulation  tests  can 
range  from  $10  to  $120  for  1-D  calculations  and  from  $100  to  $1000  for  2-D 
calculations  although  more  complex  2-D  problems  can  be  encountered  which 
may  cost  many  times  this  amount.  These  relatively  low  costs  would  indicate 
the  practicality  of  increased  computer  calculations  and  modeling  for  the 
design,  prediction,  and  analysis  of  nuclear  airblast  simulation  experiments 
which  use  high  explosives. 

Recommendations  for  future  work  include  the  development  of  equations  of  state 
and  "burn"  routines  for  the  nonideal  explosive  materials  involved  in  simulator 
explosion  processes.  Specific  calculations  of  well-defined  and  controlled 
experiments  will  be  required  to  validate  these  improvements.  The  continuation 
of  the  development  of  the  1  1 /2-D  version  of  SAP  should  be  given  high  priority 
so  calculations  of  problems  having  nearly  1-D  geometry  can  be  calculated 
without  long  and  costly  2-D  computations.  Finally,  calculations  of  specific 
explosive  configurations  currently  being  used  in  the  ongoing  research  in 
nuclear  blast  simulator  technology  should  be  accomplished  using  the  appropriate 
calculational  tools  to  provide  guidance  to  designers  of  simulation  experiments. 
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APPENDIX  A 

COMPUTER  LISTING  OF  ONE-DIMENSIONAL  SIMILARITY  SOLUTION 
PROGRAM  FOR  A  PLANAR  BLAST  WAVE  PROBLEM 
(BASIC  LANGUAGE) 
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APPENDIX  B 

HYDROCODE  CALCULATION  DATA 
REPORTS  PRODUCED  DURING  THIS 
TECHNICAL  EFFORT 


TABLE  B-l.  CONTINUED. 


906.1031 

906.1041  AST-30  SAP— Lagrangian  HAVE  HOST  T-series  Incident  wave 


TALBE  B-l .  CONCLUDED. 


APPENDIX  C 
PROGRAM  SH0K2B 
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OUTLINE  FOR  PROGRAM  SH0K2B 
SHOCK  TUBE  INTERIOR  BALLISTICS 


Note:  Region  subscripts  refer  to  Figure  56. 
Given:  Pi,  P4t  pi,  p4,  ei,  e4,  Ti,  T4,  and  TIME 


1 .  Solve  for  y i »  Y2»  ai ,  a4 


where 


Y 


P_ 

pe 


+  1 


2.  Iterative  solution*  for  P2 


.  P4  P2 

where  ^ 


1  - 


(Y4  -  1) 


&■)  k  - 


V2y7  •  2Y1  +  (y i  +  1)7? 2 


£  -  >)J 


-2y4 

T^TT 


3.  Complete  solution  for  states  in  Regions  2  and  3 

where  P3  =  P2 

Ui  =  U4  =  0 

U2  =  U3  = 


(Y4-1)  "I 

r 2a^  1 

-(! 

3\  Y  4 

(Y4  -  1)  J 

7) 

P  3  =  p4 


(fr) 


Y  4 


*Iuzzolino,  H.J.,  "Subroutine  ZEROIN,"  Air  Force  Weapons  Laboratory,  Kirtland 
Air  Force  Base,  New  Mexico,  September  1970. 
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5.  Solve  for  position  of  contact  surface 


where  X  =  U,  •  TIME 
c 


6.  Solve  for  position  of  head  and  tail  of  expansion  wave 


where 


(y,  +  1) 


•  TIME 
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(NTES,  R.  Henny) 

(NTED,  Lt  Williams,  Off.  Record  Cy.) 
(NTY) 

(HO,  R.  Minge) 

HQ  Space  Div.,  P.0.  Box  92060,  Worldway 
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(XRTB) 
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(MNNH,  M.  Delveccio) 

(MNNH,  D.  Gage) 

HQ  Space  Div.,  P.0.  Box  92060,  Worldway 
Postal  Ctr.,  Los  Angeles,  CA  90009 
(DYS) 

Foreign  Technology  Div.,  Wright-Patterson 
AFB,  OH  45433 

(NICD  Library) 

AFML,  Wright-Patterson  AFB,  OH  45433 
(Tech.  Lib.) 
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AF  Geophysics  Lab.,  Hanscom  AFB,  MA  01731 
(SUOL  Rsch .  Lib.) 

AF  Tech.  Appl .  Ctr,  Patrick  AFB,  FL  32925 
(TAP) 

(B.  Webster) 

AF  Off.  of  Sci.  Rsch.,  Bolling  AFB,  DC  20332 
(L.  Ormend) 

(W.  Best) 

AF  Inst,  of  Tech.,  Wright-Patterson  AFB, 

OH  45433 

(Library) 

AU,  Maxwell  AFB,  AL  36112 
(AUL/LDE) 

Vela  Seismological  Ctr.,  DAF,  312  Montgomery 
St.,  Alexandria,  V A  22314 
(W.  Ullrich) 

USAF  Academy,  Colorado  Springs,  CO  80840 
(DFSLB) 

DCS  for  Rsch.,  Dev.,  &  Acq.,  DA,  Washington, 
DC  20310 

(Tech.  Lib.) 

Chief  of  Engrs.,  DA,  Washington,  DC  20314 
( DAE N- RDM) 

(DAEN-MCE-D) 

DCS  for  Ops.  &  Plans,  DA,  Washington,  DC 
20310 

(Tech.  Lib.) 

NATO  School  (SHAPE),  APO  New  York,  NY  09172 
(U.S.  Documents  Officer) 

Harry  Diamond  Labs.,  DA,  Adel  phi,  MD  20783 
(DELHD-NP) 

(DELHD-TI ) 

U.S.  Army  Engr.  Waterways  Exper.  Sta. 
Vicksburg,  MS  39180 
(Tech.  Lib.) 

(J.  Jackson 
(L.  Ingram) 

(W.  Flathau) 

(R.  Walker) 

U.S.  Army  Mat.  &  Mechanics  Rsch.  Ctr. 
Watertown,  MA  02172 
(Tech.  Lib) 

U.S.  Army  Engr.  Ctr.,  Ft.  Bel  voir,  VA  22060 
(ATSEN-SY-L) 


U.S.  Army  Mat.  Dev.  &  Readiness  Cmd. 

Alexandria,  VA  22333 
(Tech.  Lib.) 

U.S.  Army  Nuclear  &  Chem.  Agency 

Springfield,  VA  22150 
(Library) 

(W.  Simms.) 

U.S.  Army  Comb.  Arms  Comb.  Dev.  Acty. 

Ft.  Leavinworth,  KS  66027 
(G.  Steger) 

BMD  Advanced  Tech.  Ctr.,  HuntsvilTe,  AL  35807 
(CRDABH-S) 

(1CRDABH-X) 

Redstone  Scientific  Info.  Ctr.,  U.S.  Army 

Missile  Cmd.,  Redstone  Arsenal,  AL  35809 
(Ch.,  Documents) 

U.S.  Army  Ballistic  Rsch.  Labs.,  Aberdeen 

Proving  Ground,  MD  21005 
(E.  Baicy) 

(W.  Taylor) 

(J.  Keefer) 

U.S.  Atm y  Engr.  Div.,  Huntsville,  AL  35807 
(HNDED-SR) 

U.S.  Army  Engr,  Div.,  Ohio  River,  Cincirn&ti, 

OH  45201 

(Tech.  Lib.) 

U.S.  Army  Mobility  Equip.  R&D  Ctr. 

Ft.  Bel  voir,  VA  22060 
(Tech.  Lib.) 

U.S.  Army  Training  &  Doctrine  Cmd.,  Ft.  Monroe 

Ft.  Monroe,  VA  23651 
(K.  Auveduti) 

U.S.  Army  Mat.  Cmd.,  Dover,  NJ  07801 
(DRCPM-NUC) 

Department  of  the  Army,  Washington,  DC  20315 
(Ch.  of  Eng.  (ENGMC-EM) ) 

Chief  of  Naval  Material,  Washington,  DC  20360 
(MAT  0323) 

Chief  of  Naval  Ops.,  Washington,  DC  20350 
(OP  981) 

Chief  of  Naval  Research,  Arlington,  VA  22217 
(N.  Perrone) 

(Tech.  Lib.) 
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Battalion  Ctr.,  Port  Hueneme,  CA  93041 
(S.  Takahashi) 

(Tech.  Lib.) 

David  W.  Taylor  Naval  Ship  R&D  Ctr. 
Bethesda,  MD  20084 
(Library) 

Naval  Facilities  Engr.  Cmd. 

Washington,  DC  20390 
(Tech.  Lib.) 

(Code  04R) 

(Code  03A) 

Naval  Postgraduate  School,  Monterey,  CA 
93904 

(Tech.  Lib.) 

Naval  Research  Laboratory,  Washington,  DC 
2037b 

(Tech.  Lib.) 

Naval  Ship  Engr.  Ctr.,  Washington,  DC  20362 
(Tech.  Lib.) 

Naval  Surface  Weapons  Center 
Silver  Spring,  MD  20910 
(Navy  Nuc.  Prgms.) 

Naval  Surface  Weapons  Center 
Dahlgren,  VA  22448 
(Tech.  Lib.) 

Naval  Weapons  Evaluation  Facility 
Kirtland  AFB ,  NM  87117 
(Tech.  Lib.) 

NW0  (Code  753),  China  Lake,  CA  93557 
(Tech.  Lib.) 

Strategic  Systems  Project  Office,  Dept,  of 
the  Navy,  Washington,  DC  20376 
(Tech.  Lib.) 

Development  Center,  Fire  Support  Branch 
Dept,  of  the  Navy,  Quantico,  VA  22134 
(L.  Gapenski) 

Naval  Ship  R&D  Center,  Bethesda,  MD  20084 
(Library 

Naval  Ship  R&D  Ctr.,  Underwater  Explosive 
Rsch.  Div.,  Portsmouth,  VA  23709 
(Tech.  Lib.) 


Office  of  Naval  Rsch.  Branch  Office 
Pasadena,  CA  91101 
(E.  Weinburg) 

Department  of  Energy,  Div.  of  Hq.  Services 
Lib.  Br.  G-043,  Washington,  DC  20545 
(Class.  Tech.  Lib.) 

Department  of  Energy,  Nevada  Operations 
Office,  Las  Vegas,  NV  89114 
(Tech.  Lib.) 

Department  of  Energy,  Albuquerque  Ops.  Off. 
Albuquerque,  NM  87115 

(Class.  Tech.  Lib.) 

Oak  Ridge  National  Laboratory,  Oak  Ridge, 

TN  37830 

(Tech.  Lib.) 

Sandia  National  Laboratories,  Kirtland  AFB, 
NM  87185 

(Info.  Dist.  Div.) 

(L.  Vortman) 

Sandia  National  Laboratories,  Livermore 
Laboratory,  Livermore,  CA  94550 
(Tech.  Lib.) 

Lawrence  Livermore  National  Laboratory 
Livermore,  CA  94550 

(Tech.  Info.  Dept.  L-3) 

(M.  Nordyke) 

Los  Alamos  National  Scientific  Laboratory 
Los  Alamos,  NM  87545 
(Reports  Lib.) 

(G.  Stillman) 

Central  Intel  1 igence  Agency 
Washington,  DC  20505 
(J.  Ingley) 

Dept,  of  the  Interior,  Bureau  of  Mines 
Denver,  CO  80225 
(Tech.  Lib.) 

(L.  Obert) 

U.S.  Geological  Survey,  Washington,  DC  20244 
(E.  Chao) 

U.S.  Geological  Survey,  Menlo  Park,  CA  94025 
(E.  Jackson) 

U.S.  Geological  Survey,  Flagstaff,  AZ  86001 
(D.  Roddy) 


157 


Aerospace  Corp.,  Los  Angeles,  CA  90009 
(Tech.  Info.  Services) 

Agbabian  Associates,  El  Segundo,  CA  90245 
(M.  Agbabian) 

Applied  Theory,  Inc.,  Los  Angeles,  CA 
90024 

(J.  Trulio) 

BDM  Corp.,  Albuquerque,  NM  87119 
(J.  Leach) 

BDM  Corp.,  McLean,  VA  22101 
(Tech.  Lib.) 

Boeing  Co.,  Seattle,  WA  98124 
(Aerospace  Lib.) 

(B.  Dyrdahl) 

California  Research  and  Tech.,  Inc. 
Woodland  Hills,  CA  91364 
(Tech.  Lib.) 

(K.  Kreyenhagen) 

Civil  Systems,  Inc.,  Albuquerque,  NM 
87110 

(J.  Bratton) 

Civil  Systems,  Inc.,  Midland,  TX  79702 
(S.  Melzer) 

University  of  Denver,  Colorado  Seminary 
Denver,  CO  80210 

(J.  Wisotski) 

General  American  Trans.  Corp.,  Niles,  IL 
60648 

(M.  Balcerzak) 

GE-TEMP0,  Albuquerque,  NM  87110 
(L.  Kennedy) 

GE-TEMP0,  Center  for  Advanced  Studies 
Santa  Barbara,  CA  93102 
(J.  Shoutens) 

(DASIAC) 

Higgins,  Auld  and  Assoc.,  Albuquerque,  NM 
87112 

(N.  Higgins) 

J.  H.  Wiggins  Co.,  Redondo  Beach,  CA  90277 
(J.  Collins) 


1 1 T  Research  Institute,  Chicago,  IL  60616 
(Tech.  Lib.) 

Kaman  AviDyne,  Burlington,  MA  01803 
(F.  Criscione) 

Kaman  Sci.  Corp.,  Colorado  Springs,  CO  80933 
(Lib.) 

Karagozian  and  Case,  Los  Angeles,  CA  90042 
(Lib. ) 

Lockheed  Missiles  and  Space  Co.,  Inc. 
Sunnyvale,  CA  94088 
(Tech.  Lib.) 

Lockheed  Missiles  and  Space  Co.,  Inc. 

Palo  Alto,  CA  94304 
(T.  Geers) 

McDonnell  Douglas,  Huntington  Beach,  CA 
92647 

(R.  Halprin) 

Merritt  CASES,  Inc.,  Redlands,  CA  92373 
(J.  Merritt) 

Mission  Research  Corp.,  Santa  Barbara,  CA 
93101 

(C.  Longmire) 

Nathan  M.  Newmark,  Consult.  Engrg.  Services 
University  of  Illinois,  Urbana,  IL  61801 
(N.  Newmark) 

Pacific-Sierra  Research  Corp.,  Santa  Monica, 
Ca  90404 

(H.  Brode) 

Pacifica  Technology,  Del  Mar,  CA  92014 
(R.  Allen) 

Physics  International  Co.,  San  Leandro,  CA 
94577 

(Tech.  Lib.) 

(F.  Sauer) 

Rand  Corp,  Santa  Monica,  CA  90406 
(Tech.  Lib.) 

R  &  D  Associates,  Marina  Del  Rey,  CA  90291 
(R.  Port) 

(J.  Lewis) 

(J.  Carpenter) 
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Science  Applications,  Inc.,  La  Jolla,  CA 
92038 

(Tech.  Lib.) 

Science  Applications,  Inc.,  Me  Lean,  VA 
22102 

(J.  Cockyane) 

Southwest  Research  Institute,  San  Antonio, 
TX  78284 

(W.  Baker) 

SRI  International,  Menlo  Park,  CA  94025 
(Tech.  Lib.) 

(B.  Gas  ten) 

Systems,  Science  and  Software,  Inc. 

La  Jolla,  CA  92038 
(T.  Cherry) 

(Tech.  Lib.) 

Teledyne  Brown  Engineering  (SETAC) 
Huntsville,  AL  35807 
(M.  Patel) 

Terra  Tek,  Inc.,  Salt  Lake  City,  UT  84108 
(Tech.  Lib.) 

Tetra  Tech.,  Inc.,  Pasadena,  CA  91107 
(Tech.  Lib.) 

University  of  Illinois,  Champaign,  IL  61820 
(S.  Hendron) 


TRW  Defense  and  Space  Sys.  Group 

Redondo  Beach,  CA  90278 
(B.  Sussholtz) 

(P.  Dai) 

(N.  Lipner) 

(Tech.  Info.  Center) 

TRW  Defense  and  Space  Systems  Group,  San 

Bernardino  Ops.,  San  Bernardino,  CA  92402 
(B.  Fay) 

(G.  Hulcher) 

University  of  Oklahoma,  Dept,  of  Info,  and 

Computing  Science,  Norman,  OK  73069 
(J.  Thompson) 

Eric  H.  Wang,  Civil  Engineering  Rsch.  Fac. 

Albuquerque,  NM  87131 
(N.  Baum) 

(J.  Berglund) 

(G.  Leigh) 

(D.  Calhoun) 

Weidlinger  Assoc.  Consulting  Engrs. 

New  York,  NY  10022 
(M.  Baron) 

Weidlinger  Assoc.  Consulting  Engrs. 
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